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Abstract 


$ 

A  study  of  2-D  acoustic  mode  coupling  in  range  dependent  shallow  water  with  a  trans¬ 
versely  isotropic  (TI)  bottom  is  presented.  Most  marine  sediments  exhibit  a  fair  de¬ 
gree  of  transverse  isotropy  that  can  have  a  significant  effect  on  the  signal  properties 
of  strongly  bottom  interacting  sound.  Locally,  transverse  isotropy  has  the  greatest 
effect  on  the  fundamental  and  near  fundamental  modal  overtones.  The  local  shallow 
water  TI  modes  have  reduced  amplitude  in  the  sediment  relative  to  the  corresponding 
shallow  water  modes  for  an  isotropic  bottom.  Calculations  of  mode-mode  coupling  co¬ 
efficients  for  a  range  dependent  medium  indicate  that  mode  coupling  is  more  strongly 
confined  to  modal  nearest  neighbors  for  a  TI  medium  characterized  predominantly  by 
shear  wave  anisotropy,  when  compared  to  the  corresponding  isotropic  medium.  As  the 
frequency  increases,  the  strongest  coupling  occurs  between  higher  overtones  and  also 
becomes  more  strongly  peaked  around  nearest  neighbors.  The  coupled  mode  theory  of 
Maupin[Geophys.  J.  93,  173-185  (1988)]  is  employed  to  model  the  coupling.  This  the¬ 
ory  can  treat  smooth  gradients  and  sloping  layer  boundaries  for  all  five  of  the  bottom 
elastic  moduli  in  a  TI  medium,  the  densities,  and  the  range  dependence  of  the  water 
column  itself.  This  coupled  mode  formulation  also  properly  accounts  for  the  continuity 
of  stress  and  displacement  boundary  conditions  in  an  exact  way  at  irregular  interfaces. 
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Introduction 


Acoustic  propagation  in  shallow  water  is  generally  characterized  by  strong  range  de¬ 
pendence  and  interaction  with  the  bottom.  The  range  dependence  is  both  geometric  and 
material  in  nature.  Fluctuations  in  water  depth  and  layer  thicknesses  of  bottom  sediments 
impose  geometric  variations  of  the  medium  with  range.  Range  dependent  gradients  in  elastic 
parameters  and  density  of  individual  layers  further  complicate  the  task  of  characterizing  and 
understanding  shallow  water  acoustic  propagation.  Adequate  modeling  and  understanding 
of  shallow  water  signals  that  interact  strongly  with  the  bottom  require  a  proper  treatment 
of  the  marine  sediment  and  basement  properties.  Marine  sediments  typically  exhibit  finite 
shear  wave  speeds  that  are  much  less  than  the  sound  speed  in  the  water  column^’^.  The 
vertical  gradients  of  the  shear  speed  can  be  quite  large,  and  velocity  anisotropy  is  an  almost 
universal  feature  of  marine  sediments,  with  transverse  isotropy  (TI)  being  the  most  common 
type  of  anisotropy^"®. 

Strong  range  dependence  causes  energy  in  an  initially  unidirectional  propagating  signal 
to  be  redistributed  among  forward  and  backward  discrete  and  continuum  modes.  Inclusion 
of  finite  shear  speed  in  the  sediment  is  necessary  to  model  the  Stoneley  wave  propagating  at 
the  water  sediment  interface  and  to  properly  account  for  the  component  of  transmission  loss 
due  to  conversion  to  shear  waves.  Acoustic  energy  can  be  scattered  from  the  water  column 
into  the  Stoneley  wave  by  range  dependence  of  the  water-sediment  interface^.  Assuming  that 
the  bottom  properties  are  isotropic  when  they  are  really  transversely  isotropic  can  lead  to 
underestimating  sediment  sound  speed  gradients,  and  overestimating  sediment  thickness  and 
shear  velocity®.  Also,  as  will  be  seen  below,  incorrectly  assuming  isotropy  has  a  significant 
effect  on  the  redistribution  of  acoustic  energy  through  range  dependence  induced  mode 
coupling. 

The  original  formulation  of  the  coupled  mode  equations  for  sound  propagation  in  a  range 
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dependent  ocean  effected  a  local  separation  of  the  Helmholtz  equation  for  the  velocity  poten¬ 
tial  or,  equivalently,  the  pressure  (Pierce®,  Milder^®).  The  dependent  variable  is  represented 
as  the  superposition  of  a  set  of  range  varying  basis  functions,  the  “local  modes,”  with  range 
dependent  amplitude  coefficients.  The  elements  of  this  local  basis  are  chosen  to  be  the  modes 
of  the  plane  layered  structure  that  corresponds  locally  in  terms  of  material  properties  and 
layer  thicknesses  to  the  range  dependent  structure.  This  approach  leads  to  another  second 
order  differential  equation  that  must  be  solved  to  obtain  the  range  dependent  modal  ampli¬ 
tude  coefficients.  The  right  hand  side  of  this  equation  consists  of  source  terms  quantifying 
the  strength  of  the  coupling  between  different  local  modes.  The  formulations  presented  by 
Pierce®  and  Milder^®  yield  first  and  second  order  coupling  coefficients  that  depend  on  the 
first  and  second  order  derivatives,  respectively,  with  respect  to  the  range  coordinate  of  the 
local  mode  functions. 

The  second  order  coupling  coefficients  are  cumbersome  to  deal  with  analytically.  They 
can  be  shown  to  depend  on  the  second  derivatives  and  the  squares  of  the  first  derivatives 
with  respect  to  the  range  coordinate  of  the  boundary  slopes  ^and  material  parameters.  Con¬ 
sequently,  the  second  order  coupling  coefficients  have  been  routinely  ignored  (Chwieroth  et 
al.^^,  Rutherford  and  Hawker^^,  McDanieP^,  HalP'^). 

It  is  interesting  to  note,  however,  that  the  presence  of  the  second  order  coupling  coeffi¬ 
cients  is  an  artifact  of  the  formulation.  It  is  a  consequence  of  working  with  the  Helmholtz 
equation  (a  second  order  differential  equation)  rather  than  directly  with  the  coupled  first 
order  equations  for  the  pressure  and  velocity.  Odom^^’^®  has  derived  a  local  coupled  mode 
theory  directly  from  the  field  quantities  pressure  and  velocity  that  contains  all  the  modal 
interaction  physics  in  a  single  coupling  coefficient.  This  formulation  exhibits  explicit  de¬ 
pendence  on  geometric  and  material  gradients,  and  is  mathematically  and  numerically  more 
efficient.  Maupin^^  extended  the  results  of  Refs. (15, 16)  to  take  elastic  effects  including 
anisotropy  into  account.  We  have  applied  Maupin’s  extensions  to  the  case  of  fiuid-elastic 
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media  in  order  to  examine  the  mode  coupling  in  a  realistic  shallow  water  model. 


Theory 


This  section  contains  a  fairly  brief,  self-contained  summary  of  the  coupled  mode  theory  for 
layered  fluid-elastic  media  as  developed  by  Maupin^^.  We  include  this  because  these  tech¬ 
niques  seem  to  be  generally  unfamiliar  to  the  ocean  acoustics  community.  A  particularly 
important  point  is  the  treatment  of  the  boundary  conditions  at  the  interface  between  two 
geometrically  irregular  layers.  Rutherford  and  Hawker^^  derived  corrections  to  the  eigen¬ 
functions  for  a  plane  layered  medium  that  satisfy  the  boundary  conditions  at  irregular  layer 
interfaces  to  first  order  in  the  interface  slope.  It  is,  however,  possible  to  satisfy  the  bound¬ 
ary  conditions  at  the  irregular  interface  exactly  by  transforming  inhomogeneous  boundary 
conditions  to  homogeneous  boundary  conditions  and  adding  an  additional  source  term  to 
the  governing  system  of  diff’erential  equations^^.  This  has  recently  been  rediscovered  by 
Fawcett^®  for  fluid  media.  Gillette^^  introduced  a  local  coordinate  transformation  that  leads 
to  a  solution  that  exactly  satisfies  the  boundary  conditions  for  the  case  of  a  single  perfectly 
rigid  range  dependent  boundary.  Gillette’s  problem  can  also  be  solved  with  the  local  mode 
theory  described  here  without  transforming  to  a  special  local  coordinate  system.  In  fact  the 
following  treatment  of  the  coupled  mode  problem  leads  to  a  solution  which  exactly  satis¬ 
fies  the  range  dependent  boundary  conditions  at  all  interfaces  with  no  approximations  or 
neglected  terms.  As  will  be  seen  below,  it  is  not  necessary  to  construct  depth  functions 
satisfying  boundary  conditions  involving  normal  derivatives  on  a  range  dependent  bound¬ 
ary.  This  exact  solution  is  also  numerically  tractible,  and  may  be  computed  using  any  good 
normal  mode  code  as  the  core  program.  The  treatment  is  valid  for  solid-solid  as  well  as  fluid- 
solid  and  fluid-fluid  boundaries.  It  is  also  valid  for  general  anisotropic  media.  Our  specific 
examples  are  carried  out  for  transversely  isotropic  media  with  a  vertical  axis  of  symmetry. 
First  we  treat  the  mode  coupling  due  to  range  dependence  within  the  bottom  material.  In 
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the  second  sub-section  we  review  the  additional  terms  in  the  mode  coupling  matrix  due  to 
the  range  dependence  of  the  water  column  itself  and  the  range  dependence  of  the  fluid-solid 
boundary  at  the  water-bottom  interface. 

A.  Mode  Coupling  in  Solid  Elastic  Media 

We  use  a  Cartesian  coordinate  system  in  which  the  x  axis  (  or  Xi  axis)  is  the  direction 
of  the  range  dependence,  y-axis  (or  X2  axis)  is  the  axis  along  which  there  is  no  variation 
in  medium  properties,  and  2:-axis  (or  0:3  axis)  is  the  depth  axis  and  taken  to  be  positive 
downward.  The  particle  displacement  vector  is  w  =  {wx,  Wy,  w^). 

For  the  elastic  moduli,  the  matrix  notation  of  Woodhouse^®  is  employed  such  that 

~  (1) 

Note  that  this  is  not  the  same  as  the  widely  used  abbreviated  subscript  notation  for  the 
elasticity  tensor  as  described  by  Auld^^  for  example.  The  individual  Cj/s  in  Woodhouse’s 
notation  are  matrices,  not  individual  matrix  elements  as  in  the  abbreviated  subscript  nota¬ 
tion. 

The  equations  of  motion  for  an  elastic  medium  are 

P-^  =  Viti  +  f,  (2) 

where  p  is  the  density,  f  is  an  external  force  and  the  traction  vector  t  is  defined  by 


where  w  =  {wx,'Wy,Wz)  and  ti  =  {jix,  ny,  Ti^). 

The  displacement,  Fourier  transformed  with  respect  to  y  and  t,  is  represented  as 

/+00  r+00 

/  w(x,  y,  z,  t)exp(i(py  -  ut))dydt 

-00  J  —00 


(3) 


(4) 
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where  p  is  the  spatial  wave  number  in  the  y-direction.  The  equations  of  motion  can  then  be 
written  as 


2  ^^3  j. 


(S) 


where 

dw  ,  dw  ^  ^ 

U  -  Cii  -  ipCi2W  +  Ci3-^.  (6) 

The  same  symbol  w  is  used  for  both  the  transformed  and  untransformed  displacement.  The 
subsequent  development  is  carried  out  completely  in  the  {x,  z,  p,  u)  domain,  and  there  should 
be  no  confusion.  The  use  of  i  as  both  an  index  and  as  \/— 1  should  also  be  clear  from  context. 

Introducing  the  6-component  displacement-stress  vector  u  =  (w,  where  w  is  as  de¬ 
fined  above  and  t  =  ti  =  (t^x,  Txy,  Txz),  the  equations  of  motion  can  be  written  as  the 
first  order  system  where  the  derivatives  with  respect  to  the  propagation  direction,  that  is, 
the  direction  of  the  range  dependence  of  the  structure,  appear  only  on  the  left  hand  side  of 
equation 


(7) 


The  differential  operator 


—CilCi3-i-  +  CiiCuip 


C 


-1 

11 


^  =  I  1  (8) 

~  '§z  (^33^)  +  +  §^{Q32W)  +P^Q22  “ -|-  ipC2iC{i  , 


does  not  depend  on  the  horizontal  derivatives.  The  matrices  Qij  are  defined  as 


Qij  —  Cij 


(9) 
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The  boundary  conditions  require  the  continuity  of  traction  and  displacement  across  inter¬ 
faces.  The  free  surface  condition  for  an  elastic  (fluid)  medium  is  that  the  traction  (pressure) 
vanishes  and  a  radiation  condition  is  assumed  as  2:  ^  oo.  The  interfaces  of  the  range  depen¬ 
dent  medium  are  taken  to  be  of  the  form  2:  =  h{x)  with  normal  n.  Thus  the  slope  of  m-th 
interface  can  be  written  as 


dzn 


dhr, 


—  hr, 


(10) 


dx  dx 

By  introducing  the  inclination  angle  6m  =  tan“^(ATO),  the  normal  vector  can  be  expressed  as 


n  =  sin  9i  —  cos  0k. 


(11) 


The  continuity  of  traction  T  =  ijiij  across  m-th  interface  can  be  written  as 

[T]  =  [tisin  dm  -  taCOS  9m]m 

=  -r  '.  [t3  -  hmt]m  =  0.  (12) 

yi  +  ^m 

The  square  brackets  [•],„  indicate  the  jump  of  the  enclosed  quantity  across  the  m*^  interface, 
taken  from  bottom  to  top.  Consequently,  we  have 


[tslm  =  hm[t]m-  (13) 

The  continuity  of  traction  normal  to  a  sloping  interface  is  then  equivalent  to  a  jump  in 
the  traction  along  the  vertical  axis.  The  equations  of  motion  (7)  along  with  the  interface 
boundary  conditions  (13)  and  the  free  surface  and  radiation  conditions  are  an  exact  formal 
representation  of  the  equations  for  the  displacement-stress  field  in  a  range  dependent  layered 
elastic  medium.  What  makes  a  solution  of  the  problem  difficult  is  the  inhomogeneous  form 
of  the  interface  boundary  condition  Eq.(13).  Historically  this  inhomogeneous  boundary 
condition  has  been  dealt  with  by  ignoring  the  inhomogeneity^^’^^’^®’^^,  and  replacing  it  with 
the  approximate  homogeneous  condition.  That  is,  the  condition  that  the  normal  component 
of  the  traction  be  continuous  across  interfaces  was  replaced  by  the  condition  that  only  the 
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vertical  component  of  the  traction  be  continuous.  It  was  pointed  out  by  Maupin^^  that  the 
traction  discontinuity  in  the  interface  boundary  conditions  can  be  converted  to  a  localized 
volume  force  located  along  the  interface.  This  follows  from  a  representation  theorem  for 
elastic  media  investigated  by  Burridge  and  KnopolP^.  The  resulting  equivalent  volume  force 
becomes  a  source  term  on  the  right  hand  side  of  Eq.(7),  and  the  interface  boundary  conditions 
become  homogeneous.  Equations  (7)  and  (13)  can  now  be  written  in  the  absence  of  body 
forces  as 


[t]m5{z  -  hmix)) 


with  the  interface  conditions 


(14) 


[talm  =  [w];„  =  0.  (15) 

Equations  (14)  and  (15)  are  a  very  important  result.  This  first  order  system  of  inhomo¬ 
geneous  equations  with  homogeneous  boundary  conditions  formally  describes  the  evolution 

j 

of  the  displacement-stress  fields  along  the  range  direction.  The  solution  to  this  system  will 
now  be  expressed  in  terms  of  coupled  local  modes.  These  local  modes,  defined  below,  are  the 
eigenfunctions  of  the  range  independent  medium  that  locally  share  the  same  depth  depen¬ 
dence  as  the  range  dependent  medium.  This  means  that  locally  at  some  point  xq  in  range, 
the  density  p(xo,  z)  and  the  elastic  moduli  Cij{xo,z)  are  taken  to  be  functions  of  depth  only 
so  that 


p{xo,z)=p{z)  and  Cij{xo,  z)  =  Cij{z).  (16) 

The  wave  propagation  problem  for  a  2-dimensional  range  dependent  medium  can  be 
solved  exactly  in  terms  of  the  eigenfunctions  of  the  range  independent  medium.  These 
eigenfunctions  are  the  homogeneous  solutions  to  Eq.(14)  with  homogeneous  boundary  con- 
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ditions  Eq.(15).  The  boundary  conditions  at  the  irregular  interfaces  are  satisfied  exactly  by 
including  the  effective  source  term  in  Eq.(14).  No  approximations  have  been  made. 

Local  homogeneous  solutions  of  the  equations  of  motion  (Eq.(14))  which  depend  parametri¬ 
cally  on  X  are  represented  as 


u(j:o,2:)exp{-zA;’'(a:o)3;}  (17) 

with  u  satisfying 

-  z/c’'(a:o)u’'(a:o,  z)  —  y4u'‘(a;o,  z)  (18) 

and  the  homogeneous  boundary  conditions  [w’’]^  =  0  and  [tgjm  =  0  across  interfaces.  The 
horizontal  wave  number  in  the  rr-direction  is  ^''(rro),  and  taken  to  be  real. 

The  final  definition  required  is  the  following  scalar  product  between  two  local  eigenfunc¬ 
tions  of  index  r  and  q. 


fOO 

{u!’,  u")  =  M  (w’^f  -  (19) 

Jo 

where  *  indicates  complex  conjugation.  The  scalar  product  (19)  is  Hermitian,  i.e. 

{L9)  =  {9jr.  (20) 

The  local  modes  at  fixed  values  of  frequency  and  p  are  orthogonal  with  respect  to  this  scalar 
product.  The  local  modes  are  normalized  such  that 


{n^u^}  =6rq.  (21) 

Thus,  they  all  carry  the  same  energy  flux  across  planes  x  =  constant. 

The  basic  idea  of  the  coupled  local  mode  technique  is  to  seek  a  solution  for  the  equations 
of  motion  as  a  coupled  set  of  local  modes  whose  amplitudes  and  phases  vary  with  laterally 
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varying  structure.  The  evolution  of  the  range  dependent  amplitude  determines  how  energy 
is  exchanged  between  modes  as  a  signal  propagates  through  the  medium.  The  solution  of 
the  equations  of  motion  for  the  displacement-stress  field  in  the  range  dependent  medium  is 
represented  as  the  sum  over  local  modes 

^  =  ^(V(a:)exp(-i  /  k^{C)dC)  ^  (22) 

t  }  r  \  Jo  J  [  J 


where  k{()  is  the  local  horizontal  wavenumber.  The  local  modes  satisfy  the  homogeneous 
boundary  conditions,  Eq.(15),  of  a  plane  layered  medium,  and  can  therefore  be  computed 
with  any  appropriate  normal  mode  code. 

The  derivation  of  the  evolution  equation  for  the  range  dependent  amplitude  coefficients 
Cr{x)  proceeds  in  the  same  fashion  as  previous  coupled  mode  developments.  The  represen¬ 
tation  Eq.(22)  is  substituted  into  the  equations  of  motion  Eq.(14).  The  scalar  product  of 
the  resulting  expression  is  formed  with  the  displacement-stress  vector  of  the  ^th  mode  u’, 
yielding: 


c 

ax 


(23) 


with  the  coupling  matrix 


_ _ _  ♦  / 

Bqr  =  {-(u^  -^)  +  i  /imW^*[t’’]m}exp  {k^  -  ¥  )dC, 

In  the  case  of  very  weak  range  dependence,  we  can  set 


(24) 


Bqr  =  0, 


(25) 


indicating  that  individual  modes  propagate  independently  without  interacting.  This  is  the 
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adiabatic  approximation^’ The  validity  of  the  adiabatic  approximation  requires  that 
the  medium  properties  change  very  slowly  with  range.  The  total  change  can  actually  be 
quite  large,  but  the  rate  of  change  with  range  must  be  small  enough  so  that  modes  do 
not  exchange  energy  among  each  other  and  are  able  to  adjust  their  shapes  to  the  local 
environment.  Backscattering  is  automatically  excluded  by  the  adiabatic  approximation,  as 
is  radiation  to  the  continuum.  A  medium  characterized  by  layer  thicknesses  or  material 
properties  that  change  significantly  over  a  mode  equivalent  ray  cycle  distance  will  not  be 
well  modeled  by  the  adiabatic  approximation. 

The  form  of  the  coupling  matrix  in  Eq.(24)  is  not  well  suited  to  numerical  computation 
because  of  the  presence  of  the  range  derivatives  of  the  local  mode  functions.  The  coupling 
matrix  Bqr  can  be  transformed  so  that  the  only  range  derivatives  appearing  in  the  expression 
are  of  the  density  and  elastic  moduli  within  layers  and  of  the  interface  functions  hm{x)  at 
layer  boundaries.  We  omit  the  lengthy  derivation  and  refer  the  interested  reader  to  Reference 
(17).  The  final  form  for  the  coupling  matrix  in  the  solid  medium  comprising  the  bottom  is 
then: 


Bqr  — 


A:*? 


^(r< 


g*  ■  2  T  dw'>*  ■  .  ■  dw^ 

pwV  -  -  W®  ipQ23-^ 


„  dwi* 

+  -^Q32ipW^  -  w'^*g22P  w'’  -  -  w'?*Zp(C'2lCf/)t'' 

o  T 


+  w’*g22wV  +  )t'’ 


exp  {i  [(k’ -k’')dC; 


(26) 


The  coupling  matrix  Eq.(26)  involves  an  integral  term  related  to  the  lateral  derivative 
of  the  elastic  moduli  and  density  inside  the  layers.  The  interface  term  is  a  combination  of 
an  expression  derived  from  the  continuity  conditions  and  another  arising  from  jumps  in  the 
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lateral  derivatives  of  the  elastic  coefficients.  It  is  the  former  of  the  these  two  interface  terms 


S<V 


that  appears  as  the  effective  volume  source  term  on  Eq.(14).  There  are  no  range  derivatives 
of  the  local  eigenfunctions,  and  there  is  no  need  to  introduce  a  special  coordinate  system. 
The  expression(26)  describes  the  coupling  in  a  fully  anisotropic  2-dimensional  medium.  We 
have  derived  the  form  of  the  matrices  Cij  and  Qij  for  a  transversely  isotropic  medium,  and 
listed  the  result  in  the  Appendix.  These  have  also  been  given  by  Maupin^®.  We  use  the 
notation  of  Love^®  for  the  five  elastic  moduli  required  to  characterize  a  transversely  isotropic 
medium.  The  relationship  between  Love’s^®  notation  and  the  abbreviated  subscript  notation 
for  the  elastic  moduli  of  a  TI  medium  is  also  given  in  the  Appendix. 

B.  Coupling  Terms  for  Fluid  Layer  and  Fluid-Elastic  Interface 

In  this  subsection  we  very  briefly  review  the  the  additional  terms  required  to  treat  the 
mode  coupling  due  to  the  range  dependence  of  the  water  column  and  the  fluid-solid  interface 
at  the  water-bottom  boundary^^.  We  make  the  two  assumptions  that:  1.  The  water  column 
surface  is  flat;  and  2.  The  water  density  p  and  incompressibility  A  may  vary  with  depth  and 
range  but  not  with  time.  We  assume  an  ocean  of  depth  h  —  h{x). 

In  a  fluid  medium,  the  displacement-stress  vector  u  is 


The  X  component  of  the  particle  displacement  is  lOa,,  and  the  scalar  t  is  the  pressure.  The 
equations  of  motion  written  in  first  order  matrix  form  are 


^=Au  +  f  (28) 

where  the  operator  A  in  the  fluid  is 
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0 


(29) 


i_-2L  +  A 
A  pa;2  '  dz 


(_L  A) 

\  pLo"^  dz  ) 


.  -poj  0 

and  f  is  the  source  vector.  As  before  p  is  the  spatial  wavenumber  in  the  y  (cross-range) 
direction. 

An  ideal  inviscid  fluid  does  not  support  shear,  so  the  boundary  condition  at  the  fluid- 
solid  bottom  boundary  requiring  continuity  of  the  tangential  component  of  displacement  is 
relaxed  and  replaced  with  a  free  slip  boundary  condition.  This  free  slip  boundary  condition 
results  in  a  physically  unrealistic  discontinuity  in  the  tangential  component  of  displacement 
Wx,  which  can  be  remedied  by  including  a  small  nonzero  viscosity  in  the  equations  of  motion 
for  the  fluid. 

The  continuity  conditions  across  a  fluid-solid  interface  of  the  form  z  —  h{x)  are 


[ta  —  Ati]  =  0,  \wz  —  hwx]  =  0. 


(30) 


Again,  these  continuity  conditions  can  be  transformed  into  continuity  conditions  across 

j 

a  flat  interface  plus  a  volume  force  located  along  the  fluid-solid  interface.  A  first  part  of  this 
force  applies  in  the  solid  part  only: 


f  =  /i[-ti]5(2;  —  h{x)). 

A  second  one  can  be  taken  as  applying  on  either  side  of  the  interface; 


(31) 


fi  =  h  1  -  K^)) 

where  has  to  be  taken  in  the  layer  where  we  actually  choose  to  apply  the  force. 
The  scalar  product  (19)  now  becomes; 


(32) 


rh  foo 

(u^  u'-)  =  i  /  {wl*f  -  f^*wl)dz  +  i  /  -  tTw^)dz 

Jh 


(33) 
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To  obtain  an  evolution  equation  for  the  amplitude  coefficients,  one  uses  the  properties  of 
the  scalar  product  (33)  and  proceeds  exactly  as  in  the  solid-solid  case.  The  transformation 
of  the  coupling  matrix  Bqr  also  proceeds  in  an  identical  manner.  The  final  result  requires 
that  several  terms  be  added  to  the  right  hand  side  of  Eq.(26).  On  the  solid  side  of  interface 
it  is  necessary  to  add 


-■A(*:«-r)«r;.+r«>a  (34) 

to  the  terms  which  already  hold  for  a  purely  solid  medium.  On  the  fluid  side  of  the  interface 
it  is  necessary  to  add  the  term 


P 


f  + 


(35) 


\  dz  dz  ~  \X  puj^ 

Finally,  inside  the  fluid  layer  are  the  integral  terms  describing  the  coupling  resulting  from 
the  continuous  range  dependent  variations  of  the  water  column 


i  df  ^ 

~  r  +  - -—dz. 


/o  \A  p(^^  J  dz  pu^  dz 

These  additional  terms  are  multiplied  by  the  global  term: 


(36) 


^;^^exp(i£(i«-fcOrfc)  (37) 

A  final  note  is  that  the  integral  in  Eq.(26),  describing  coupling  due  to  continuous  range 
dependence  of  the  density  and  elastic  moduli  in  the  bottom,  now  extends  from  the  ocean 
bottom  h{x)  to  infinity.  It  would  also  be  a  simple  matter  to  incorporate  an  additional  solid 
layer  at  the  surface  to  model  the  mode  coupling  in  an  ice  covered  sea. 
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C.  Energy  Conservation 


As  indicated  by  Eq.(21),  the  local  modes  are  normalized  to  carry  the  same  energy  flux 
across  planes  x  —  constant.  We  can  obtain  a  statement  of  energy  conservation  for  a  lossless 
range  dependent  medium  by  substituting  the  local  mode  representation  of  the  displacement- 
stress  field  Eq.(22)  into  the  scalar  product  Eq.(19)  and  setting  the  derivative  with  respect 
to  X  equal  to  0,  yielding 

=  E  +  B;,)  Cr(rt)c;(l) 

Q,r 

=  0.  (38) 

We  have  used  the  fact  that 

j 

E  (E  B*rCr{x)^  =  E  (E  ^r9C*(x)  j  ,  (39) 

since  the  mode  indices  {q,  r}  are  summed  over  the  same  set. 

The  only  way  for  Eq.(38)  to  be  satisfied  generally  is  for  the  coupling  matrix  to  be  anti- 
Hermitian,  i.e. 


=  -b;,.  («) 

This  anti-Hermiticity  is  a  necessary  consequence  of  energy  conservation  in  a  lossless  medium. 
It  can  be  seen  by  inspection  that  the  coupling  matrix  B^r  given  by  Eq.(26)  is  anti-Hermitian. 
The  additional  coupling  terms  for  the  fluid  layer  and  interface,  equations  (34),  (35)  and  (36) 
are  also  anti-Hermitian. 
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An  obvious  consequence  of  Eq.(40)  is  that 


w 


Re{Brr)  -  0.  (41) 

If  we  insist  that  the  phase  of  the  local  modes  be  continuous  from  point  to  point  in  the 
medium  then  we  should  also  choose 

Im{Brr)  -  0.  (42) 

D.  Evolution  of  Range  Dependent  Mode  Amplitudes 

The  theory  we  have  summarized  up  to  this  point  only  describes  the  local  mode-mode  in¬ 
teractions  caused  by  the  range  dependence  of  the  medium.  In  order  to  synthesize  a  complete 
signal  propagating  in  a  strongly  range  dependent  medium  for  which  the  adiabatic  approxi¬ 
mation  is  not  valid,  we  must  solve  the  evolution  equation  Eq.(23)  for  the  mode  amplitudes. 
A  complete  solution  for  a  strongly  range  dependent  medium  must  account  for  the  interaction 
of  both  forward(-l-)  and  backward(-)  propagating  modes.  The  evolution  of  both  the  forward 
and  backward  propagating  components  of  a  signal  is  described  by 

Ic-wj 

where  c"*"  and  c“  are  n-dimensional  vectors  whose  elements  are  the  amplitude  coefficients  of 
n  forward  and  backward  propagating  modes.  The  nxn  matrices  B++,  B“+  and  B 

describe  forward-to-forward,  forward-to-backward,  backward-to-forward  and  backward-to- 
backward  coupling,  respectively. 

If  we  specify  a  geometry  defined  by  a  heterogeneous  region  sandwiched  between  two 
homogeneous  (plane  layered)  regions,  and  assume  a  signal  incident  from  the  left  onto  the 
heterogeneous  region,  Eq.(43)  defines  a  2nx2n  boundary  value  problem  for  the  amplitudes  of 


'B''”''(a:)  B+  (x)\  /  C^lx) 
,B— '■(a;)  B  {x)  j  \c“(a;), 


(43) 
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the  forward  and  backward  propagating  modes.  A  schematic  of  a  range  dependent  medium  is 
shown  in  Figure  1.  The  boundary  values  are  c^{xl),  known  at  a;  =  on  the  left  side  of  the 
heterogeneous  region,  and  c~{xr)  =  0  on  the  right  side  of  the  heterogeneous  region  dXx  —  xr. 
Stable  numerical  solution  of  the  two  point  boundary  value  problem  defined  by  Eq.(43)  and  the 
two  boundary  values  is  problematic  due  to  presence  of  growing  and  decaying  mode  amplitudes 
within  the  heterogeneous  region.  This  situation  is  exacerbated  if  the  heterogeneous  region 
is  extended  in  range.  Kennett^^  reformulated  the  two  point  boundary  value  problem  as 
an  initial  value  problem  for  the  refiection  and  transmission  properties  of  the  heterogeneous 
region  using  invariant  imbedding  techniques. 

The  procedure  is  to  define  a  transmission  matrix  T{xr,xl)  that  connects  the  c'^{xl)  on 
the  left  side  of  the  heterogeneous  region  with  the  on  the  right  side  of  the  heterogeneous 

region 


=  T(xr,  xr)c+{xl).  (44) 

In  addition  a  reflection  matrix  R(x/e,  xr)  is  defined  that  relates  the  backscattered  component 
c~(xi)  from  the  heterogeneous  region  to  the  forward  propagating  component  c'^(xl)  at  the 
left  side  of  the  heterogeneous  region 

c~(xl)  =  R(xr,Xr)c^(xl).  (45) 

We  differentiate  Eqs.(44)  and  (45)  with  respect  to  xr 

(47) 
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The  derivatives  of  the  amplitude  vectors  are  replaced  with  their  expressions  from  Eq.(43),  and 
c'~{xl)  can  be  removed  from  the  equation  using  Eq.(45).  After  removing  a  common  factor 
of  we  arrive  at  coupled  matrix  Ricatti  equations  for  the  reflection  and  transmission 

matrices  for  the  heterogeneous  region 


and 


dxr 


R  =  B-+  +  B— R  -  RB++  -  RB+-R 


(48) 


The  integration  begins  to  the 
left  with  the  initial  conditions 


_d_ 

dxi 


T  =  -TB++  -  TB+~T. 


(49) 


right  of  the  heterogeneous  region  at  Xr  and  proceeds  to  the 


T(a;/j,XH)  =  I  and  R(xij,  x;?)  =  0.  (50) 

I  is  the  n  X  n  identity  matrix. 

j 

The  stability  of  the  coupled  set  of  first  order  equations  (48)  and  (49)  is  significantly 
improved  over  that  of  the  original  two  point  boundary  value  problem.  If  backscattering  from 
the  heterogeneous  region  can  be  neglected,  then  we  are  left  with  only  the  equation 

^T=-TB-.  (51) 

Eq.(51)  describes  a  level  of  approximation  between  that  implied  by  the  adiabatic  approxi¬ 
mation  Eq.(25)  and  the  complete  solution  described  by  Eq.(48)  and  Eq.(49)  with  the  initial 
conditions  Eq.(50). 

This  concludes  our  summary  of  the  complete  2-D  coupled  mode  theory.  This  theory 
incorporates  an  exact  treatment  of  the  boundary  conditions  at  sloping  interfaces  as  well  as 
being  able  to  handle  smooth  gradients  in  density  and  all  elastic  moduli.  We  have  applied 
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the  theory  to  a  transversely  isotropic  medium  with  a  vertical  symmetry  axis.  The  theory 
is,  however,  valid  for  general  anisotropic  media  at  the  expense  of  additional  algebraic  and 
numerical  complexity.  In  a  general  anisotropic  medium  the  horizontally  polarized  shear 
waves  do  not  decouple  from  the  compressional  waves  and  vertically  polarized  shear  waves. 
Consequently,  the  motion  is  represented  by  a  single  6x6  system  of  equations  rather  than  the 
separate  4  x  4(P-SV)  and  2  x  2(SH)  systems  required  for  isotropic  and  transversely  isotropic 
media. 


Numerical  Results 

In  this  section  we  numerically  investigate  the  effects  of  transverse  isotropy  on  the  modal 
dispersion,  eigenfunctions  and  the  coupling  matrix  B^r  for  a  realistic  shallow  water  model. 
We  are  in  the  process  of  coding  the  matrix  Ricatti  equations  (48)  and  (49),  and  will  present 
results  from  the  integration  of  those  equations  in  a  future  paper. 

Elastic  anisotropy  is  a  well  established  geoacoustic  property  of  marine  sediments^"®.  The 
most  common  form  of  anisotropy  observed  in  marine  sediments  is  transverse  isotropy(TI) 
with  a  vertical  axis  of  symmetry.  In  a  TI  medium  propagation  in  all  azimuthal  directions 
lying  in  planes  including  the  symmetry  axis  is  equivalent.  In  TI  marine  sediments  with  ver¬ 
tical  symmetry  axes,  horizontally  polarized  and  horizontally  propagating  waves  travel  faster 
than  vertically  propagating  and  vertically  polarized  waves.  The  primary  mechanisms  for 
transverse  isotropy  in  marine  sediments  are:  1.  aligned  cracks  and  pores;  2.  recrystallization 
of  anisotropic  minerals;  and  3.  compositional  layering  on  a  very  fine  scale.  If  isotropy  is 
erroneously  assumed  sound  speed  gradients  are  underestimated,  sediment  layer  thicknesses 
can  be  overestimated  and  shear  velocity  can  be  overestimated®.  The  most  likely  mechanism 
for  the  observed  transverse  isotropy  of  sediments  at  the  water-bottom  interface  is  composi¬ 
tional  layering.  Recrystallization  of  anisotropic  minerals  may  be  important  for  deeper  lying 
sediments,  but  to  depths  of  at  least  lA:m  beneath  the  ocean  bottom,  compositional  layering 
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is  believed  to  be  the  dominant  mechanism®. 

We  have  examined  the  effect  of  transverse  isotropy  on  the  mode  structure  and  mode 
coupling  of  a  bottom  interacting  signal  propagating  in  shallow  water.  The  code  used  to 
generate  the  eigenvalues,  eigenfunctions  and  kinetic  energy  integrals  for  a  plane  layered 
fluid-elastic  TI  medium  was  DISPERSO^®.  The  computations  were  carried  out  on  a  Sun 
SPARCstation  LX,  mostly  in  double  precision.  At  the  highest  frequencies  it  was  necessary 
to  employ  quadruple  precision  for  the  first  few  eigenvalues  and  eigenfunctions.  Our  shallow 
water  model  was  taken  from  Berge  et  al.^^  who  analyzed  multi-component  seismic  data  from 
an  experiment  conducted  in  21m  of  water  about  10km  east  of  New  Jersey. 

The  data  analyzed  by  Berge  et  al.^®  were  collected  in  1986  by  Roundout  Associates  Inc. 
and  Woods  Hole  Oceanographic  Institution.  Initial  attempts  at  modeling  the  data  assuming 
isotropy  of  the  bottom  material  were  not  considered  adequate^®,  and  led  to  further  modeling 
efforts  employing  an  anisotropic  reflectivity  program.  The  resulting  TI  models  estimated 
by  Berge  et  al.^®  provided  a  good  fit  to  their  data  and  constrained  four  of  the  five  elastic 
parameters  necessary  to  describe  a  TI  medium.  They  did  not  have  enough  resolution  to 
determine  the  compressional  wave  anisotropy  and  assumed  Cn  =  C33  (A  =  C  in  Love’s^® 
notation). 

In  order  to  specify  the  anisotropy,  we  introduce  the  dimensionless  parameters 


and 


V 


=  1-77 


1  - 


A-2L’ 


=  1-j. 

Values  of 


(52) 


(53) 
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T)'  —  0  and  —  0 


(54) 


indicate  isotropy.  The  departure  from  isotropy  increases  as  rj'  and/or  (j)'  change  from  0. 

The  dynamic  stability  of  a  transversely  isotropic  medium  requires  that  certain  conditions 
be  imposed  on  the  elastic  moduli  (Auld^^,  Postma^°,  Backus^^): 


A>N>0, 

C>  0, 

L  >  0, 

F‘^<C{A-N).  (55) 

The  most  plausible  physical  mechanism  for  TI  with  a  vertical  symmetry  axis  within  the 
interfacial  sediments  is  fine  compostional  layering.  This  mechanism  imposes  the  additional 
constraints 


C>F, 

c>Il, 

N>L  (56) 

(Backus^^).  A  final  independent  inequality  constraint  for  compositionally  layered  TI  media 
is 


{A~L){C-L)>{F  +  L)^  (57) 

(Postma^°,  Berryman^^). 

The  model,  shown  in  Figure  2,  consists  of  a  21m  thick  isovelocity  water  layer  over  12.5m 
of  TI  sediments,  followed  by  12.5m  of  isotropic  sediments  with  a  steep  gradient  in  both 
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compressional  and  shear  speeds.  The  density  of  the  sediments  is  taken  to  be  constant  at 
2100A:^/m^.  The  TI  layer  is  characterized  by  r]'  —  0.012  and  (j)'  =  0.  The  base  of  the  model 
is  an  isotropic  elastic  half-space  with  a  shear  speed  of  1450m/s,  a  compressional  speed  of 
3000m/s,  and  a  density  of  2A00kg/m^. 

We  have  explicitly  focused  on  the  qP-qSV  component  of  propagation  in  the  sediments. 
The  excitation  of  horizontally  polarized  shear  waves(SH)  in  the  sediments  by  a  propagating 
acoustic  signal  in  the  water  requires  3-D  heterogeneity  or  anisotropy  of  a  more  general  nature 
than  we  treat. 

Figure  3  illustrates  the  phase  velocity  dispersion  for  the  model  of  Figure  2.  The  modal 
phase  velocities  for  the  TI  model  are  generally  higher  than  for  the  corresponding  isotropic 
medium.  This  was  apparently  first  noted  by  Stoneley^^,  and  is  a  consequence  of  greater 
material  stiffness  sampled  by  components  of  wave  particle  motion  parallel  to  the  bedding 
plane.  The  difference  between  the  modal  phase  velocities  for  the  TI  and  isotropic  media 
increases  with  increasing  frequency.  As  the  frequency  increases,  the  phase  velocity  of  the  TI 
mode  can  approach  the  phase  velocity  of  the  next  higher  isotropic  mode.  This  is  also  evident 
in  the  group  velocity  curves  shown  in  Figure  4. 

Despite  a  relatively  small  difference  in  the  fundamental  mode  phase  velocity  between 
the  TI  and  isotropic  models,  the  eigenfunctions  are  substantially  different.  Because  the 
eigenfunction  calculation  depends  on  the  inverse  of  the  eigen-phase  velocity,  small  changes 
in  the  smallest  eigen-phase  velocity  can  have  significant  effects  on  the  computation  of  the 
corresponding  eigenfunction.  This  is  also  the  reason  for  using  quadruple  precision  at  the 
higher  frequencies  for  the  first  few  modes.  DISPER80  directly  integrates  the  equations  of 
motion  which  become  quite  stiff  at  high  frequencies.  Efficient  methods  of  generating  the 
fluid-elastic  modes  are  important  for  high  frequency  applications. 

The  eigenfunctions  for  the  vertical  component  of  displacement  at  lOHz  and  20 Hz  are 
shown  in  Figures  5  and  6,  respectively.  Although  it  is  not  visible  on  the  scale  at  which  the 
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modes  are  plotted,  all  eigenfunctions  have  been  normalized  such  that  the  vertical  component 
of  displacement  is  unity  at  the  sea  surface.  The  scale  of  the  other  three  components  of  the 
eigenfunctions  is  derived  from  this  normalization.  This  choice  of  normalization  is  somewhat 
arbitrary  and  merits  discussion  since  the  energy  integrals  and  their  partial  derivatives  (Fig¬ 
ures  7(a-c),  8  and  9)  are  computed  using  this  normalization.  The  purpose  of  Figures  5  and  6 
is  to  illustrate  the  effect  of  the  transverse  isotropy  on  the  eigenfunctions.  The  eigenfunctions 
could  have  been  normalized  so  the  isotropic  and  corresponding  TI  eigenfunctions  had  equal 
energy,  or  so  that  they  had  equal  peak  amplitudes  or  equal  amplitudes  at  the  water  sediment 
interface.  Although  we  have  not  computed  the  Green’s  function  for  the  medium,  which  forces 
the  selection  of  an  explicit  source  location,  we  are  assuming  that  the  source  and  receiver  are 
located  within  the  water  column  or  at,  or  slightly  beneath,  the  water-sediment  interface.  Our 
chosen  normalization  emphasizes  differences  near  the  water-sediment  interface.  Modes  with 
a  given  surface  amplitude  will  have  very  different  amplitudes  at  the  water-sediment  interface 
depending  on  whether  the  bottom  is  isotropic  or  anisotropic.  Had  we  chosen  a  normalization 
such  that  the  peak  vertical  displacement  was  unity,  corresponding  amplitudes  in  the  water 
column  would  be  quite  different.  Our  normalization  yields  eigenfunctions  with  different  en¬ 
ergies  in  isotropic  and  corresponding  TI  media,  but  because  they  have  the  same  vertical 
displacement  at  the  surface,  the  differences  in  energy  reflect  the  different  amount  of  energy 
input  into  the  medium  required  to  produce  the  same  surface  displacement.  The  contribution 
of  an  individual  mode  to  the  Green’s  function  for  the  medium  is  directly  proportional  to 
the  mode  amplitude  at  the  chosen  source  depth.  Therefore  our  choice  of  normalization  will 
emphasize  differences  resulting  from  shallow  sources. 

Comparison  of  eigenfunctions  between  isotropic  and  anisotropic  media  is  somewhat  prob¬ 
lematic  in  any  case,  since  there  is  no  simple  method  of  establishing  correspondence  between 
the  two.  Because  of  the  constraint  relations  among  the  elastic  parameters  listed  earlier,  it  is 
not  possible  to  arbitrarily  perturb  one  elastic  modulus  without  corresponding  perturbations 
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to  other  moduli.  Arbitrary  perturbations  to  elastic  moduli  can  destroy  the  symmetry  of  the 
elastic  stiffness  tensor  and  produce  unphysical  results  in  calculations. 

In  the  lQ-2QHz  frequency  range  the  water  depth  is  approximately  A/4,  which  means  that 
the  bottom  is  near  an  acoustic  radiation  maximum.  A  propagating  acoustic  signal  in  this 
band  will  thus  be  dominated  by  the  fundamental  mode  guided  along  the  water  sediment 
interface.  For  frequencies  greater  than  approximately  IHz  for  the  isotropic  medium  and 
llHz  for  the  TI  medium,  the  phase  velocity  of  the  fundamental  mode  is  less  than  the 
sediment  shear  speed  at  the  water  sediment  interface.  The  fundamental  mode  above  the 
threshhold  frequency  is  therefore  a  Stoneley  wave.  At  frequencies  lower  than  the  threshhold, 
the  mode  is  more  properly  termed  a  pseudo-Rayleigh  wave. 

It  can  be  clearly  seen  from  Figures  5  and  6  that  the  fundamental  carries  the  most  energy  in 
this  frequency  band.  The  generally  smaller  peak  amplitude  of  the  low  order  TI  modes  relative 
to  the  isotropic  modes  is  a  consequence  of  the  greater  material  stiffness  for  horizontally 
propagating  waves  in  the  TI  medium.  As  the  mode  number  increases  for  a  given  frequency, 
the  eigenfunctions  persist  to  greater  depths  in  the  structure.  The  fractional  amount  of  modal 
energy  within  the  10m  thick  TI  layer  decreases,  so  the  influence  of  the  anisotropy  on  the 
eigenfunctions  also  decreases.  Also,  as  the  frequency  of  a  given  mode  increases  the  component 
of  the  mode  in  the  bottom  becomes  more  and  more  like  pure  SV.  It  can  be  seen  from  Figure 
3  that  as  the  frequency  of  a  mode  increases,  its  phase  velocity  approaches  the  shear  speed 
at  the  sediment  interface.  The  compressional  speed  is  much  higher  than  the  sediment  shear 
speed,  and  so  the  compressional  component  of  the  mode  is  evanescent,  leading  to  the  almost 
pure  SV  behavior  of  the  modes  at  high  frequency. 

Figures  7(a-c)  shows  the  relative  mode  kinetic  energy  as  a  function  of  frequency  for  the 
first  three  modes.  The  kinetic  energy  Em  of  the  m***  mode  is 
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(58) 


The  energy  in  the  Stoneley  wave  peaks  at  around  llHz.  Berge  et  al.^®  have  plotted  amplitude 
vs.  frequency  for  the  unfiltered  vertical  component  of  the  data  from  one  of  their  profiles. 
(Their  Figure  7,  and  here  reprodnced  as  our  Figure  7d.)  The  main  features  of  the  lower 
frequency  part  of  their  spectrum  are  well  represented  by  our  kinetic  energy  plots.  The 
Stoneley  wave  peak  at  about  llHz  in  Berge  et  al.’s^®  data  appears  to  be  particularly  well 
modeled  by  our  mode  calculations.  The  low  amplitude  maxima  appearing  at  1.8Hz  and  bHz 
in  Figures  7a  and  7b,  respectively,  and  the  kink  at  approximately  8Hz  in  Figure  7c  occur 
at  the  knees  of  the  modal  dispersion  curves(Figure  3).  These  knees  occur  at  the  frequency 
at  which  a  mode  becomes  trapped  in  the  sediment  layer.  The  mode  dies  out  very  rapidly  in 
the  underlying  halfspace,  and  the  phase  velocity  drops  abruptly  towards  the  sediment  shear 
wave  speed.  The  peak  at  A.5Hz  in  Figure  7c  corresponds  to  the  very  sharp  minimum  in  the 
group  velocity  curve  for  m  =  2  (Figure  4).  At  this  point  the  group  velocity  rises  to  meet 
the  phase  velocity  at  the  cutoff  value.  The  comparisons  between  our  model  calculations  and 
Berge  et  al.’s^®  data  is  not  direct  as  we  have  not  attempted  to  correct  our  modeling  results 
in  such  a  way  that  would  permit  absolute  amplitude  comparisons  with  their  experimental 
data. 

Berge  et  al.^®  reported  that  their  synthetic  seismograms  were  quite  sensitive  to  small 
pertnrbations  in  the  elastic  modulus  F  (ci3  in  the  abbreviated  subscript  notation).  The 
modulus  F  affects  the  the  propagation  of  qP  and  qSV  at  angles  intermediate  between  the 
horizontal  and  vertical  in  a  TI  medium  with  a  vertical  symmetry  axis.  In  an  attempt  to 
illuminate  the  sensitivity  of  a  bottom  interacting  acoustic  signal  to  sediment  anisotropy,  we 
computed  partial  derivatives  of  mode  energy  with  respect  to  a  parameter  r]' . 

The  dimensionless  partial  derivative  of  the  mode  energy  is  defined  by 
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1  dEm  _  I  dF  dEjn  _  A -2L  dEm 
'K'W  ~  ^  Em  dF  ■ 

In  Figure  8  the  partial  derivative  Eq(59)  is  plotted  versus  t]'  x  100  for  the  first  three  modes 
at  lOHz.  Since  a  value  of  =  0  indicates  isotropy,  the  anisotropy  increases  with  increasing 
values  of  the  abcissa.  The  magnitudes  of  the  derivatives  for  the  first  three  modes  are  relatively 
large  and  negative.  The  negativity  of  the  derivatives  indicates  that  the  energy  of  the  modes 
will  decrease  with  increasing  anisotropy  as  measured  by  increasing  rj'.  Again,  we  mention  that 
the  vertical  component  of  the  eigenfunctions  are  all  normalized  to  have  unit  displacement  at 
the  surface. 

Physically  Figure  8  illustrates  that,  as  the  elastic  modulus  F  departs  from  its  isotropic 
value,  less  energy  is  required  to  produce  the  same  vertical  surface  displacement.  (In  an 
isotropic  medium,  F  corresponds  to  the  Lame  parameter  A.)  Our  perturbation  reduces  F 
from  its  isotropic  value,  thereby  reducing  the  stiffness  of  the  medium  somewhat  at  angles  in¬ 
termediate  between  the  horizontal  and  vertical.  The  effect  is  the  same  as  reducing  the  spring 
constant  of  a  mass-spring  system.  Less  energy  is  required  to  produce  a  given  displacement  of 
the  mass  suspended  from  a  weaker  spring.  The  relatively  large  dimensionless  magnitudes  of 
the  derivatives  are  a  measure  of  the  sensitivity  of  the  eigenfunction,  and  hence  the  acoustic 
signal,  to  anisotropic  medium  perturbations.  Since  the  fundamental  has  the  largest  deriva¬ 
tive,  it  will  also  be  most  sensitive  to  changes  in  F.  Although  not  shown  in  Figure  8,  at  lOHz 
for  the  mode  m  =  3,  dEm/ dr]'  is  a  very  weak  function  of  rj'  and  nearly  zero.  The  mode 
m  =  3  persists  to  greater  depth  into  the  bottom,  and  has  relatively  more  energy  both  in  the 
water  column  and  the  underlying  halfspace. 

We  also  found  for  this  model  that  a  2.4%  change  in  77'  could  produce  a  15%  change  in  the 
phase  velocity  of  the  Stoneley  wave  mode.  Over  the  range  of  r]'  from  0  to  0.024,  the  phase 
velocity  of  the  fundamental  increases  from  145.57m/s  to  169.67m/s. 
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The  frequency  sensitivity  of  the  mode  energy  for  a  single  value  of  the  anisotropy  parameter 
r]'  =  0.012  is  illustrated  in  Figure  9  by  a  plot  of  dEm/drj'  versus  frequency.  The  main  thing 
to  notice  is,  that  over  the  4Hz  to  2()Hz  frequency  range  depicted,  the  magnitude  of  the 
derivative  is  increasing.  A  shallow  water  signal  with  a  strongly  excited  Stoneley  wave  will 
become  more  sensitive  to  the  anisotropy  with  increasing  frequency.  Of  course  as  the  frequency 
increases  still  further,  the  excitation  of  the  Stoneley  wave  will  drop  off,  and  the  influence  of 
the  anisotropy  on  that  part  of  the  signal  will  also  decrease. 

We  have  not  computed  the  corresponding  partial  derivatives  with  respect  to  (j)',  which 
controls  the  P  wave  anisotropy,  while  holding  r)'  =  0.  This  combination  of  anisotropy 
parameters  violates  the  condition  expressed  by  Eq.(57),  while  the  converse,  i.e.  0'  =  0  and 
r]'  >  0  does  not.  We  have,  however,  computed  an  example  that  shows  the  effect  of  both 
ri'  >  0  and  0'  >  0.  As  previously  mentioned,  Berge  et  al.^®  did  not  investigate  the  effects 
of  P-wave  anisotropy  because  their  profiles  were  not  long  enough  to  resolve  it.  The  length 
of  their  long  profile  was  about  200m.  The  P-wave  anisotropy  could  be  important  for  longer 
propagation  paths,  since  values  as  high  as  40%  have  been  reported  in  shale^^,  although  8% 
to  14%  is  probably  more  typical  of  marine  sediments^. 

Figure  10  is  a  plot  of  the  vertical  displacement  component  of  the  fundamental  eigenfunc¬ 
tion  for  the  isotropic  case  rj'  —  0  and  0'  =  0  (solid  line),  and  for  the  two  transversely  isotropic 
cases  corresponding  to  rj'  =  0.012  and  (/>'  =  0  (dashed  line),  and  r/'  =  0.012  and  (j)'  =  0.012 
(dotted  line).  A  value  of  77'  =  0.012  corresponds  to  Berge  et  al.’s^®  long  profile  model  (their 
Table  1).  The  vertical  components  of  the  eigenfunctions  are  all  normalized  to  have  unit 
displacement  at  the  surface.  The  amplitude  of  the  fundamental  for  the  case  rj'  —  0.012  and 
0'  =  0.012  lies  between  the  the  amplitude  for  the  isotropic  case  and  for  the  case  77'  =  0.012 
and  0'  =  0.  The  inclusion  of  this  very  modest  amount  of  P-wave  anisotropy  draws  the 
appearance  of  the  fundamental  back  towards  isotropy.  It  may  be  possible  to  increase  the 
P-wave  anisotropy  enough  so  that  there  is  essentially  no  detectable  difference  between  the 
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shapes  of  the  isotropic  and  anisotropic  eigenfunctions,  but,  in  reality  the  P-wave  anisotropy 
near  the  water-sediment  interface  is  likely  to  be  even  smaller  than  the  value  used  for  our 
calculations®.  The  phase  velocity  of  the  fundamental  for  the  case  r)'  =  0.012  and  0'  =  0  is 
159.83m/s,  which  is  actually  slightly  higher  than  the  horizontally  propagating  qSV  speed 
(158.26m/s). 

Our  ultimate  goal  is  to  improve  our  ability  to  model  and  predict  shallow  water  acoustic 
signal  propagation  in  a  range  dependent  medium.  We  have  computed  the  coupling  matrix 
Bqr  (Eq.  26)  including  the  fluid-solid  boundary  interaction  terms  described  in  Eq.(34)  and 
Eq.(35)  for  the  isotropic  model  {t}'  =  0  and  <j)'  —  0)  equivalent  to  Berge  et  al.’s  model  (Fig. 
11a, b);  for  the  TI  model  of  Berge  et  al.^^  (77'  =  0.012  and  (j)'  —  0)  (Fig.  llc,d);  and  for  a 
model  {rf  =  0.012  and  (j)'  =  0.012)  incorporating  the  weak  P-wave  anisotropy  (Fig.  lle,f). 
The  calculations  were  done  at  two  frequencies  WHz  and  20Hz. 

Berge  et  al.’s^®  TI  model  is  a  range  independent  model.  What  we  have  computed  is 
the  coupling  matrix  for  a  model  with  the  same  local  vertical  structure  as  Berge  et  al.^® 
The  absolute  values  of  the  coupling  matrix  in  Figure  TI  represent  the  effects  of  the 
geometric  medium  properties  only.  The  layer  boundary  slopes,  /i„  =  dhn/dx,  have  been  set 
equal  to  1.  All  material  parameter  horizontal  gradients  such  as  p  have  been  set  equal  to  0. 
The  absolute  values  of  the  elements  of  all  six  coupling  matrices  depicted  in  Figure  11  have 
been  normalized  by  dividing  all  elements  by  the  largest  matrix  element  of  the  entire  set  so 
comparisons  can  be  made  between  frequencies  and  between  the  isotropic  and  both  TI  media. 

As  the  frequency  increases,  the  shallow  water  structure  can  support  a  greater  number  of 
modes,  so  there  are  more  modes  to  participate  in  the  coupling.  There  is  a  preferred  mode 
pair  for  which  the  strongest  coupling  occurs  in  Berge  et  al.’s^^  model.  At  10  Hz  this  is  mode 
pair  {1, 2}  (Fig.  11c),  and  at  20  Hz  the  strongest  interaction  occurs  for  pair  {4, 5}  (Fig.  lid). 
Coupling  strength  decreases  away  from  the  diagonal  and  also  away  from  the  preferred  mode 
pair.  The  coupling  of  the  preferred  pair  is  stronger  in  the  TI  medium  than  in  the  equivalent 
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isotropic  medium.  This  latter  effect  illustrates  the  most  striking  difference  between  the  TI 
and  isotropic  media,  and  is  a  qualitative  indicator  of  the  difficulties  that  may  be  encountered 
by  ignoring  sediment  anisotropy  when  attempting  to  model  range  dependent  shallow  water 
acoustic  propagation.  Away  from  the  diagonal,  other  mode  pairs  participate  nearly  equally 
in  the  coupling  process.  This  indicates  that  a  careful  examination  of  the  coupling  matrix  is 
necessary  before  deciding  on  a  truncated  mode  set  to  employ  when  synthesizing  complete 
propagating  signals  in  a  range  dependent  medium.  Use  of  too  small  a  mode  set  will  alias 
the  coupling,  and  affect  the  amplitude  and  phase  of  a  synthesized  signal,  but  also,  it  may 
occur  that  coupling  is  confined  to  a  small  number  of  model  configurations  leading  to  more 
efficiency  in  calculations. 

We  have  also  computed  the  coupling  matrix  for  a  TI  medium  containing  weak  P-wave 
anisotropy  in  addition  to  the  S-wave  anisotropy  of  Berge  et  al.’s^®  model  (Fig.  lle,f).  The 
addition  of  weak  P-wave  anisotropy  to  the  model  actually  makes  the  medium  appear  more 
isotropic  for  the  case  shown.  It  is,  however,  important  to  remember  that  Figure  11  represents 
just  a  local  snapshot  of  coupling  interactions.  To  get  a  complete  picture,  the  actual  range 
dependent  structure  of  the  medium  must  be  imposed  and  Eqs.(48)  and  (49)  for  the  reflection 
from  and  transmission  through  the  heterogeneous  region  must  be  solved.  Because  of  the 
almost  order  of  magnitude  difference  between  the  shear  and  compressional  wavelengths  at 
the  water-sediment  interface,  the  effects  of  the  shear  wave  anisotropy  will  accumulate  much 
faster  for  a  strongly  bottom  interacting  shallow  water  acoustic  signal. 


E. Summary  and  Conclusions 

We  have  summarized  a  coupled  mode  theory  for  fluid-elastic  media  that  is  formulated  as 
a  coupled  set  of  first  order  equations,  and  accounts  exactly  for  the  inhomogeneous  boundary 
condition  due  to  range  dependent  interface  irregularities’^'^.  We  have  applied  this  theory  to 
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a  realistic  shallow  water  model  derived  from  experimental  data.  A  particular  feature  of  this 
work  is  the  inclusion  of  transversely  isotropic  bottom  sediments. 

Our  modeling  results  show  that  there  can  be  significant  qualitative  and  quantitative  dif¬ 
ferences  between  the  eigenvalues  and  modes  in  shallow  water  models  with  an  isotropic  and 
a  tranversely  isotropic  bottom.  These  differences  are  also  reflected  in  the  mode  coupling 
induced  by  range  dependence.  The  Stoneley  wave  at  the  water  sediment  interface  is  partic¬ 
ularly  sensitive  to  the  transverse  isotropy  of  the  sediments.  Conversion  to  Stoneley  waves 
has  been  shown  to  be  an  important  loss  mechanism  by  Hawker^®.  In  light  of  the  sensitivity 
of  the  Stoneley  waves  to  the  transverse  isotropy  of  the  bottom,  and  the  apparent  ease  with 
which  they  can  be  excited  by  bottom  roughness’',  some  care  should  be  taken  regarding  inter¬ 
pretations  of  strongly  bottom  interacting  acoustic  signals  derived  from  models  that  assume 
isotropic  sediment  properties. 
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Appendix:  Elastic  Moduli  Matrices  for  a  Transversely  Isotropic  Medium 


We  list  the  elastic  moduli  matrices  and  Qij  used  in  the  calculation  of  the  coupling 
matrix  Bgr,  Eq.(26). 


C'r/C'12  - 


Cfi'Cia 


/1/A  0  0  \ 

0  1/N  0 

Vo  0  1/Lj 
/O  {A-2N)/A  0\ 
1  0  0 
VO  0  0/ 

/O  0  F/A\ 

0  0  0 

Vl  0  0  / 


The  Qij,  defined  as 


Qij  —  Cij  Ci\C]^i  Cij, 


are: 


Q22  — 


/O  0  0\ 

0  iN{A-N)/A  0 

Vo  0  lJ 


/O  0 


0  \ 


Q23  =  1  0  0  2NF/A 
VO  L  0  J 
/O  0  0 


Q33  — 


0  L  0 

Vo  0  {AC-F^)/Aj 
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When  A  =  C  —  \  2/i,  L  —  N  =  fi,  and  F  =  X  the  medium  is  isotropic. 

The  relationship  between  the  notation  of  Love^®  and  the  abbreviated  subscript  notation  for 
the  elastic  moduli  of  a  TI  medium  is 


A  —  Cn,  F  —  Ci3, 


F/  C33  j  L  C44  j  TV  cqq  . 
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Figure  Captions 


Figure  1.  The  figure  illustrates  the  geometry  of  the  range  dependent  medium  modeled 
by  the  coupled  Ricatti  equations  Eq.(48)  and  Eq.(49).  The  range  dependent  region  (II)  is 
sandwiched  between  two  plane  layered  (range  independent)  regions  (I)  and  (III)  that  need 
not  be  the  same.  The  integration  of  Eq.(48)  and  Eq.(49)  proceeds  backwards  from  the  point 
xr  to  the  point  xl- 

Figure  2.  The  transversely  isotropic  shallow  water  model  of  Berge  et  al.^®  (Their 
Tablel.).  The  12.5m  thick  sediment  layer  immediately  below  the  water  sediment  inter¬ 
face  is  transversely  isotropic.  We  have  terminated  their  model  with  an  isotropic  half-space 
with  a  shear  speed  of  1450m/s,  a  compressional  wave  speed  of  3000m/s  and  a  density  of 
2Amkglm^. 

Figure  3.  Phase  velocity  dispersion  curves  for  the  shallow  water  niodel  of  Fig.  2. 
Note  that  the  phase  velocities  in  the  TI  medium  are  generally  higher  than  in  the  equivalent 
isotropic  medium.  As  the  frequency  increases,  the  phase  velocity  of  a  TI  mode  can  approach 
the  phase  velocity  of  the  next  higher  isotropic  mode.  This  is  also  true  of  the  group  velocities 
plotted  in  Fig.  4. 

Figure  4.  Group  velocities  of  the  first  three  modes  for  the  model  shown  in  Figure  2. 

Figure  5.  The  vertical  displacement  eigenfunctions  for  the  first  four  modes  at  lOHz. 
Notice  that  for  mode  3,  the  isotropic  and  TI  vertical  displacements  are  virtually  identical, 
but  that  the  amplitude  of  the  TI  displacement  in  the  halfspace  is  slightly  greater  than  the 
isotropic  displacement.  Modes  are  normalized  to  have  unit  vertical  displacement  at  the  water 
surface.  At  lOHz,  the  signal  is  dominated  by  the  fundamental  (Stoneley)  mode. 

Figure  6.  The  vertical  displacement  eigenfunctions  for  the  first  four  modes  at  20Hz.  The 
isotropic  and  TI  modes  are  distinct.  In  fact,  although  not  shown  here,  there  are  significant 
differences  between  the  isotropic  and  TI  modes  for  the  first  five  modes. 
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Figure  7.  Kinetic  energy  as  a  function  of  frequency  for  the  first  three  modes  (a-c). 
These  match  quite  well  the  lower  frequency  part  of  the  spectrum  of  some  shallow  water 
data  (d)  shown  by  Berge  et  al.^®  (Their  Figure  7,  used  by  permission  of  Blackwell  Scientific 
Publications  Ltd.)  The  Stoneley  wave  peak  at  about  11  Hz  in  Berge  et  al.^®  appears  to 
particularly  well  modeled  by  our  mode  calculations.  Note  that  the  frequency  ranges  of  (a-c) 
and  (d)  are  not  the  same. 

Figure  8.  Nondimensional  partial  derivatives  of  mode  energy  with  respect  to  the  pa¬ 
rameter  rf.  The  negativity  of  the  derivatives  indicates  that  the  mode  energy  for  the  first 
three  modes  decreases  upon  departure  from  isotropy. 

Figure  9.  The  partial  derivative  of  the  fundamental  mode  energy  with  respect  to  77' 
evaluated  at  r)'  =  0.012  as  a  function  of  frequency.  Note  that  the  absolute  value  of  the 
magnitude  of  the  derivative  generally  increases  over  the  plotted  frequency  range.  Indicating 
an  increasing  sensitivity  with  frequency  of  the  fundamental  mode  to  the  anisotropy.  A  value 
of  rf  =  0.012  corresponds  to  the  best  fitting  TI  model  for  the  long  profile  data  of  Berge  et 
al.’s^^  (their  Table  1). 

Figure  10.  Vertical  displacement  component  of  the  fundamental  eigenfunction  for  the 
isotropic  medium(solid  line,  phase  velocity=  145.57m/s),  and  for  TI  media  characterized  by 
T]'  =  0.012(dashed  line,  phase  velocity=  159.83m/s)  and  rj'  =  (j)'  —  0.012(dotted  line,  phase 
velocity  =  153.73m/s).  They  have  all  been  normalized  to  unit  displacement  at  z  =  0.  The 
frequency  is  IQHz.  The  speed  of  horizontally  propagating  qSV  waves  {y/l/p)  at  the  water 
sediment  interface  is  158.26m/s. 

Figure  11.  The  absolute  value  of  the  elements  of  the  coupling  matrix  Bgr  Eq.(26)  for 
the  model  in  Fig.  2.  Layer  boundary  slopes  h  have  been  set  equal  to  1,  and  all  horizontal 
material  parameter  gradients  such  as  p  have  been  set  to  0.  This  emphasizes  the  eflPects  of 
geometric  (boundary)  heterogeneities.  The  absolute  values  of  the  elements  of  the  six  coupling 
matrices  have  been  normalized  by  the  largest  matrix  element  of  the  entire  set  so  comparisons 
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can  be  made  between  frequencies  and  between  isotropic  and  TI  media.  The  array  elements 
have  been  normalized  so  that  dark  red  is  unit  coupling  and  dark  blue  is  zero  coupling.  The 
diagonals  have  been  purposely  left  blank  to  reflect  our  choice  of  phase  for  the  local  modes 
(Eq.(42). 
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SUMMARY 


Elastic  wave  propagation  in  strongly  scattering  media  is  described  by  an  energy  diffusion 
equation  derived  directly  from  the  coupled  mode  equations.  These  equations  provide  an 
exact  representation  of  the  displacement-stress  field  in  2-D  heterogeneous  media.  By  allow¬ 
ing  the  mode  spacing  to  approach  zero,  the  continuum  limit  of  the  modal  representation  of 
the  elastic  wavefield  is  derived.  The  final  expressions  require  neither  the  location  of  closely 
spaced  eigenvalues  nor  construction  of  eigenfunctions.  All  assumptions  and  approximations 
must  be  explicitly  stated  and  implemented  in  the  derivation  of  the  energy  diffusion  equation. 
Strong  forward  scattering  is  assumed  to  dominate.  The  energy  diffusivity  is  a  range  and  fre¬ 
quency  dependent  functional  of  the  displacement-stress  field  components  and  the  horizontal 
gradients  of  the  medium  properties  including  anisotropy.  The  diffusivity  functional  is  de¬ 
rived  directly  from  the  continuum  limit  of  the  mode  coupling  matrix;  it  is  essentially  the 
spatial  autospectrum  of  the  coupling  matrix  weighted  by  a  function  describing  the  density 
of  modes  in  spectral  space.  The  approach  presented  in  this  paper  is  in  contrast  to  energy 
diffusion  equations  derived  from  radiative  transfer  theory  in  which  the  the  diffusivity  must 
be  specified  separately  in  an  ad  hoc  manner.  The  diffusivity  functional  is  computable  from 
any  wavenumber  integration  type  code  that  generates,  or  can  be  modified  to  generate,  the 
depth  dependent  Green’s  function  for  a  fluid-elastic  medium.  Although  our  energy  diffusion 
equation  is  range  dependent,  the  computation  of  the  diffusivity  depends  on  local  medium 
properties  and  field  values  for  a  plane  layered  medium.  An  additional  difference  between 
the  diffusion  equation  of  this  paper  and  previously  published  treatments  of  elastic  energy 
diffusion  is  that  this  paper  describes  energy  diffusion  in  spectral  space,  e.g.  slowness,  as  a 
function  of  range  rather  than  diffusion  in  physical  space  as  a  function  of  time. 

Key  words:  coupled  local  modes,  energy  diffusion,  forward  scattering,  random  media, 
surface  waves. 


1 


1  INTRODUCTION 


Deterministic  models  are  quite  successful  for  modeling  the  large  and  intermediate  scales  of 
the  structure  of  the  earth.  The  scales  of  the  deterministic  problems  are  defined  by  inhomo¬ 
geneities  that  are  large  compared  to  the  wavelength  of  the  seismic  or  acoustic  waves.  Initially 
coherent  signals  tend  to  retain  their  coherence  over  long  propagation  distances.  When  the 
scale  of  the  inhomogeneities  is  small  compared  to  the  wavelength  and  are  randomly  dis¬ 
tributed,  propagating  signals  are  scattered,  and  become  progressively  more  incoherent  as 
the  range  from  the  source  increases. 

A  number  of  different  approaches  have  been  used  to  model  the  scattered  wavefield  de¬ 
pending  on  the  spatial  density  and  properties  of  the  inhomogeneities  in  the  medium.  For 
widely  spaced  scatterers,  the  Born  approximation,  which  does  not  conserve  energy,  has  been 
used  (Wu  and  Aki  1985a,  Snieder  1986).  (The  definition  of  “widely  spaced”  depends  on 
the  mean  free  path,  correlation  length,  and  frequency.)  For  more  closely  spaced  scatterers 
various  multiple  scattering  formulations  have  been  proposed  (e.g.  Wu  1985).  Phase  screen 
methods  have  recently  been  extended  to  to  the  P-SV  problem  (Fisk  &  McCartor  1991;  Fisk, 
Charette  &  McCartor  1992).  Phase  screen  methods  have  also  been  combined  with  the  Born 
approximation  to  model  propagation  in  media  with  large  numbers  of  imbedded  scatterers 
(Wu  1994).  When  the  scatterer  density  is  so  high  that  all  the  energy  is  scattered,  phase  in¬ 
formation  is  lost  and  the  field  is  completely  incoherent.  In  this  case  radiative  transfer  theory 
has  been  applied  to  derive  an  energy  diffusion  equation.  The  application  of  radiative  transfer 
theory  to  seismic  problems  was  first  suggested  by  Wesley  (1965),  followed  by  applications 
to  lunar  seismograms  by  Dainty  et  al.  (1974),  Dainty  and  Toksoz  (1977),  Nakamura  (1977) 
and  Oberst  (1985).  Wu  (1985)  applied  radiative  transfer  theory  to  study  the  coda  of  local 
earthquakes.  In  this  paper  a  new  method  of  modeling  energy  diffusion  in  strongly  scattering 
media  is  described.  Although  we  also  arrive  at  a  diffusion  equation  for  the  elastic  energy. 
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there  are  a  number  of  differences  between  the  theory  presented  in  this  paper  and  radiative 
transfer  theory.  The  main  differences  are  summarized  in  the  next  three  paragraphs. 

Radiative  transfer  theory  is  a  phenomenological  theory  that  starts  with  an  assumed 
seismic  energy  density  defined  as  an  integral  over  all  space.  Substitution  of  the  energy 
density  into  a  continuity  equation  derived  from  energy  conservation  leads  to  a  diffusion 
equation  describing  the  redistribution  of  energy  in  space  and  time.  The  derivation  of  the 
diffusion  equation  depends  on  a  number  of  implicit  as  well  as  explicit  assumptions.  The 
diffusivity,  which  should  contain  details  about  the  physical  properties  of  the  scatterers,  is 
specified  in  an  ad  hoc  manner  based  on  plausibility  arguments  about  the  size,  distribution 
and  scattering  strength  of  the  imbedded  inhomogeneities  and  partitioning  of  energy  between  . 
The  model  that  is  either  explicitly  or  implicitly  assumed  is  generally  of  a  random  distribution 
of  scatterers  imbedded  in  a  homogeneous  background.  An  advantage  of  radiative  transfer 
theory  over  the  theory  presented  in  this  paper  is  that  it  is  inherently  a  three  dimensional 
theory,  and  backscattering  is  automatically  accounted  for. 

We  derive  an  energy  diffusion  equation  directly  from  the  exact  coupled  mode  represen¬ 
tation  of  the  solution  to  the  elastic  wave  equation.  This  makes  the  theory  of  this  paper 
a  physically  based  theory  rather  than  a  phenomenological  theory.  The  assumptions  and 
approximations  must  be  explicitly  stated  in  order  to  extract  a  diffusion  equation  from  the 
starting  coupled  mode  representation.  One  key  assumption  is  that  backscattering  is  ne¬ 
glected.  Because  the  scattered  energy  is  assumed  to  be  forward  scattered,  our  theory  is  most 
appropriately  applied  to  scattering  from  velocity  heterogeneities  which  affect  the  wavefield 
phase.  Backscattering  produced  by  strong  impedance  heterogeneities  is  not  modeled  by  this 
approach  (Wu  &  Aki  1985a).  Whereas  previously  published  treatments  of  elastic  energy 
diffusion  model  energy  diffusion  in  space  and  time,  we  model  energy  diffusion  in  spectral 
space  as  a  function  of  propagation  range. 

As  shown  below,  the  diffusivity  arises  from  the  continuum  limit  of  the  mode  coupling 
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matrix,  and  it  contains  all  the  interaction  physics  and  medium  properties  so  that  the  back¬ 
ground  medium  and  imbedded  inhomogeneities  are  treated  on  an  equal  footing.  The  coupling 
between  the  compressional  and  shear  waves  is  automatically  included.  The  frequency  de¬ 
pendence  of  the  diffusion  process  and  correlations  of  the  horizontal  gradients  of  the  elastic 
moduli,  density  and  layer  boundaries  appear  explicitly  in  the  diffusivity.  All  the  required 
quantities  are  readily  computable. 

The  coupled  mode  equations  for  an  elastic  medium  are  reviewed  in  Section  2.  From 
the  coupled  mode  equations,  a  coupled  energy  equation  is  derived,  and  then  reduced  to  an 
energy  diffusion  equation.  The  derivation  is  2-D  and  valid  for  general  anisotropic  fluid-elastic 
media. 

In  Section  3.  an  approximate  solution  of  the  coupled  mode  equations  in  the  forward  scat¬ 
tering  limit  is  described,  to  be  used  later  in  the  derivation  of  the  coupled  energy  equations. 

In  Section  4.  energy  conservation  in  a  range  dependent  medium,  and  its  relation  to  the 
anti-Hermiticity  of  the  mode  coupling  matrix  is  discussed. 

In  Section  5.  the  coupled  energy  equations  are  derived  from  the  coupled  mode  equations. 
The  phase  information  carried  by  individual  modes  is  relinquished  for  a  description  of  the 
evolution  of  the  average  energy  of  a  mode  propagating  in  a  range  dependent  medium. 

In  Section  6.  the  energy  diffusion  equation  is  derived  from  the  coupled  energy  equations. 
The  energy  diffusivity  functional  is  shown  to  be  essentially  the  spatial  autospectrum  of  the 
coupling  matrix  weighted  by  the  local  density  of  modes  in  the  random  medium.  The  diffusion 
equation  can  be  transformed  to  a  Schrodinger  equation,  which  may  offer  some  advantages 
for  numerical  implementation.  However,  the  form  of  the  potential  is  somewhat  complicated. 

In  Section  7.  an  initial  condition  and  physically  reasonable  choices  for  two  boundary 
conditions  are  discussed.  The  boundary  conditions  must  be  satisfied  in  spectral  space,  rather 
than  physical  space. 

In  Section  8.  solutions  to  two  elementary  spectral  diffusion  problems  are  presented  and 
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discussed.  The  solutions  are  represented  as  Fourier  series  of  the  energy  spectrum  in  — 
space.  The  homogenization  of  the  spectrum  with  increasing  propagation  distance  through 
the  scattering  medium  is  readily  apparent  upon  examining  the  elementary  solutions.  The 
final  two  sections  include  a  brief  discussion  of  attenuation  and  our  conclusions. 


2  THE  COUPLED  MODE  EQUATIONS 

Seismic  wave  modeling  in  fluid-elastic  media  with  variations  in  2  and  3  dimensions  is  com¬ 
plicated  by  the  non-separability  of  the  spatial  independent  variables  of  the  equations  of 
motion.  However  the  normal  mode  representation  can  still  be  used  with  modification.  A 
local  separation  of  the  equations  of  motion  is  effected  by  representing  the  wavefield  compo¬ 
nents  as  the  superposition  of  a  set  of  range  varying  basis  functions,  the  “local  modes,”  with 
range  dependent  amplitude  and  phase.  The  elements  of  this  local  basis  are  the  modes  of  the 
plane  layered  structure  that  corresponds  locally  to  the  range  dependent  structure  in  terms 
of  material  properties  and  layer  thicknesses.  The  derivation  uses  the  coupled  mode  theory  of 
Maupin  (1988),  which  is  an  extension  to  elastic  media  of  the  acoustic  coupled  mode  theory 
of  Odom  (1981,  1986).  Kennett’s  (1984)  theory  uses  the  modes  of  a  reference  structure;  this 
theory  will  suffer  in  accuracy  if  the  boundaries  of  the  reference  structure  diverge  from  those 
of  the  actual  structure.  Bostock  (1992)  has  extended  Kennett’s  reference  structure  method 
to  3-D,  and  Tromp  (1994)  has  made  a  partial  extension  of  the  local  mode  method  to  3-D. 

Maupin’s  (1988)  theory  is  summarized  here,  since  it  forms  the  basis  of  the  derivation  of 
the  energy  diffusion  equation.  The  theory  is  valid  for  solid-solid  as  well  as  fluid-solid  and 
fluid-fluid  boundaries.  It  is  also  valid  for  general  anisotropic  media. 

A  Cartesian  coordinate  system  is  used  where  the  x  axis  (  or  Xi  axis)  is  the  direction  of 
the  range  dependence,  and  2:-axis  (or  0:3  axis)  is  the  depth  axis  and  taken  to  be  positive 
downward.  The  medium  is  assumed  to  be  2-D,  laterally  heterogeneous  and  composed  of  n 
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layers.  The  layers  may  have  non-planar  interfaces  hn{x),  where  the  free  surface  is  designated 
ho{x).  Material  properties,  the  density  p{x,z)  and  stiffnesses  Cijki{x,z),  may  vary  smoothly 
in  both  the  x  and  z  directions. 

For  the  elastic  moduli,  the  matrix  notation  of  Woodhouse  (1974)  is  employed  such  that 

~  ^kilj-  (1) 

The  equations  of  motion  for  an  elastic  medium  are 

=  Viti  +  f,  (2) 

where  p  is  the  density,  f  is  an  external  force  and  the  traction  vector  t  is  defined  by 


t<  —  Ci 


'  d'w 


The  displacement  is  Fourier  transformed  with  respect  to  t,  as 

/+00 

w(a;,  2:,  t)exp{—iujt)dt. 

-00 


(3) 

(4) 


The  same  symbol  w  is  used  for  both  the  transformed  and  untransformed  displacement.  The 
equations  of  motion  for  a  2-D  heterogeneous  medium,  rewritten  to  isolate  the  derivatives 
with  respect  to  x,  the  propagation  direction,  on  the  left  hand  side  are 


5u 

dx 

with  the  interface  conditions 


Au  +  J2hn 

n 


"  1 

[t]nS{z  -  hn{x))  J 


(5) 


[tsjn  =  [w]„  ==  0.  (6) 

The  square  brackets  [•]„  indicate  the  jump  of  the  enclosed  quantity  across  the  interface, 
taken  from  bottom  to  top.  The  free  surface  condition  for  an  elastic  (fluid)  medium  is  that 
the  traction  (pressure)  vanishes  and  a  radiation  condition  is  assumed  as  2  — 00.  The 
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4-component  vector  u  =  (w-jt)"^  where  w  =  {wx,Wz)  and  t  =  ti  =  {txx,  Txz),  and  the 
differential  operator  A  is 

(  -cr.’cuf  cf/  \ 

^  =  .  (7) 

V -  K  (033£)  -fCaiCfi'l 

and  Qij  —  ^ij 

Continuity  of  the  traction  normal  to  a  sloping  interface  in  the  laterally  heterogeneous 
medium  is  equivalent  to  a  jump  in  traction  along  the  vertical  axis.  Maupin  &  Kennett 
(1987)  transformed  the  traction  discontinuity  in  the  interface  boundary  conditions  to  a  lo¬ 
calized  volume  force  located  along  the  interface  by  applying  a  representation  theorem  for 
elastic  media  (Burridge  &  Knopoff  1964).  The  resulting  equivalent  volume  force  becomes  the 
source  term  on  the  right  hand  side  of  eq.  (5),  and  the  interface  boundary  conditions  become 
homogeneous. 

Eqs.  (5)  and  (6)  are  a  first  order  system  of  inhomogeneous  equations  with  homogeneous 
boundary  conditions  that  formally  describes  the  evolution  of  the  stress-displacement  fields 
along  the  propagation  direction.  The  solution  to  this  system  is  expressed  in  terms  of  coupled 
local  modes.  These  local  modes,  defined  below,  are  the  eigenfunctions  of  the  range  indepen¬ 
dent  medium  that  locally  share  the  same  depth  dependence  as  the  range  dependent  medium. 
Hence  locally  at  some  point  Xq  in  range,  the  density  p{xo,  z)  and  the  elastic  moduli  Cij{xo,  z) 
are  functions  of  depth  only  so  that 

p{xo,z)=p{z)  and  Cij{xo,z)  -  Cij{z).  (8) 

The  wave  propagation  problem  for  a  2-dimensional  range  dependent  medium  can  be 
solved  exactly  in  terms  of  the  local  eigenfunctions  of  the  range  independent  medium.  These 
eigenfunctions  are  the  homogeneous  solutions  to  eq.  (5)  with  homogeneous  boundary  con- 
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ditions  eq.  (6).  The  boundary  conditions  at  the  irregular  interfaces  are  satisfied  exactly  by 
including  the  effective  source  term  in  eq.  (5).  No  approximations  have  been  made  (Maupin 
1988). 

Local  solutions  of  the  equations  of  motion  (eq.  (5))  which  depend  parametrically  on  x 
are  represented  as 


u{xo,z)exp{-ik''{xo)x)  (9) 

satisfying 

-  ik'' (xo)u'' (xo,  z)  =  An''{xo,  z)  (10) 

and  the  homogeneous  boundary  conditions  [w'‘]„  =  0  and  [tgjn  =  0  across  interfaces.  The 
horizontal  wave  number  in  the  x-direction  is  and  taken  to  be  real. 

Maupin  (1988)  introduced  the  following  Hermitian  scalar  product  between  two  local 
eigenfunctions  of  index  r  and  q. 

rOO 

(u-?,  n^)  =  i  (w''*!’-  -  t'^V’')d^  (11) 

«/ 0 

where  *  indicates  complex  conjugation.  The  local  modes  are  orthogonal  with  respect  to  this 
scalar  product  at  fixed  values  of  frequency  and  p.  The  local  modes  are  normalized  such  that 

(u'',u’')  =  (5^,.  (12) 

Thus,  they  all  carry  the  same  energy  flux  across  planes  x  =  constant. 

The  coupled  local  mode  technique  seeks  a  solution  for  the  equations  of  motion  as  a  coupled 
set  of  local  modes  whose  amplitudes  and  phases  vary  with  laterally  varying  structure.  The 
evolution  of  the  range  dependent  amplitude  determines  how  energy  is  exchanged  between 
modes  as  a  signal  propagates  through  the  medium.  The  solution  of  the  equations  of  motion 
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for  the  stress-displacement  field  in  the  range  dependent  medium  is  thus  represented  as  the 
sum  over  local  modes 

_  /  /.X  X  r  z)  ] 

=  ^Cg{x)exp(-i  k'^{^)dn  I  (13) 

9  \  Jo  /  (  t‘^(^x,z)  J 

where  k{^)  is  the  local  horizontal  wavenumber  at  a;  The  local  modes  satisfy  the 

homogeneous  boundary  conditions,  eq.  (6),  of  a  plane  layered  medium,  and  are  therefore 
readily  computed. 

The  derivation  of  the  evolution  equation  for  the  range  dependent  amplitude  coefficients 
Cr(x)  proceeds  by  substituting  the  representation  eq.  (13)  into  the  equations  of  motion  eq.  (5). 
The  scalar  product  of  the  resulting  expression  is  formed  with  the  displacement-stress  vector 
of  the  Q'th  mode  u^,  yielding: 


^  =  B,rCr,  (14) 

which  describes  the  evolution  of  the  range  dependent  amplitude  coefficients. 

The  explicit  form  for  the  coupling  matrix  Bgr  has  been  reproduced  for  reference  from 
Maupin  (1988),  and  is  listed  in  the  Appendix  (eq.  (Al)).  The  coupling  matrix  Bgr  is  com¬ 
posed  of  products  of  components  of  the  displacement-stress  field,  their  vertical  derivatives 
and  lateral  gradients  of  the  medium  properties.  There  are  no  range  derivatives  of  the  lo¬ 
cal  eigenfunctions.  The  explicit  expression  for  the  coupling  matrix  eq.  (Al)  describes  the 
coupling  in  a  fully  anisotropic  2-D  medium. 
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3  APPROXIMATE  SOLUTION  OF  THE  COUPLED  MODE  EQUATIONS 


In  this  section  we  derive  an  approximate  solution  of  the  coupled  mode  equations  that  will 
be  required  for  the  derivation  of  the  coupled  energy  equation  in  the  next  section. 

For  a  complete  description  of  propagation  in  a  general  heterogeneous  medium,  the  evolu¬ 
tion  eq.  (14)  for  the  mode  amplitudes  must  be  solved  for  both  forward(-l-)  and  backward(— ) 
propagating  modes. 


/B++(a;)  B+“(2:) 
c^{x)  j  \B“‘'‘(2:)  B  (x) 

where  c+  and  c~  are  n-dimensional  vectors  whose  components  are  the  amplitude  coefficients 
of  n  forward  and  backward  propagating  modes.  If  we  specify  a  geometry  defined  by  a 
heterogeneous  region  sandwiched  between  two  homogeneous  (plane  layered)  regions  and 
assume  a  signal  incident  from  the  left  onto  the  heterogeneous  region,  eq.  (15)  defines  a  2n  x  2n 
boundary  value  problem  for  the  amplitudes  of  the  forward  and  backward  propagating  modes. 
The  boundary  values  are  C^^xl)  known  B,tx  —  xl  on  the  left  side  of  the  heterogeneous  region, 
and  c~{xr)  =  0  on  the  right  side  of  the  heterogeneous  region  at  x  =  Xr.  Figure  1.  illustrates 
the  medium  geometry  and  various  length  scales  employed  in  this  paper. 

We  now  assume  that  energy  is  scattered  predominantly  in  the  forward  direction  so  that 
backscattering  can  be  neglected.  The  forward  scattering  assumption  also  implies  that  this 
theory  is  applicable  to  scattering  from  medium  heterogeneities  that  are  predominantly  ve¬ 
locity  heterogeneities  rather  than  impedance  heterogeneities.  Backscattering  is  controlled 
primarily  by  impedance  contrast  (Wu  &  Aki  1985a,b). 

Applying  the  forward  scattering  assumption  eq.  (15)  can  be  written 

^  =  Bp'{x)c+{x)  (16) 
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which  can  be  approximately  integrated  (Marquering  fe  Snieder  1994) 


c+  (a:)  =  exp  ^  c+  {xl).  (17) 

The  validity  of  this  approximate  solution  to  eq.  (16)  requires  that  the  commutator  of 
and  be  small.  A  derivation  of  the  approximate  solution  eq.  (17),  and  an  excellent 

discussion  of  the  subtleties  involved  is  given  by  Marquering  &  Snieder  (1994).  We  represent 
the  matrix  exponential  in  eq.  (17)  by  its  power  series  expansion  (e.g.  Golub  &  Van  Loan 
1983,  p.  396) 


OO 

exp{Agr)  =  E 

j=0  J- 

so  that  eq.  (17)  becomes 

N  Y  ^  ^ 

(^)  ~  C.q(Xi^  +  AqrCj-{xi,)  +  —  Aqj-  ArsCs{Xi)  +  .  .  . 
r=0  r=0  l-O 

where 


(18) 


(19) 


A,,=  r  (20) 

Jxl 

Marquering  &:  Snieder  (1994)  have  used  the  same  power  series  representation  of  the  matrix 
exponential  (eq.  (19))  to  derive  a  very  efficient  algorithm  for  the  solution  of  the  coupled  mode 
equations  in  the  forward  scattering  limit  that  takes  second  order  coupling  interactions  into 
account.  In  our  derivation  of  the  coupled  energy  equations  below,  we  retain  only  the  first 
two  terms  of  eq.  (19). 


11 


4  ENERGY  CONSERVATION 


As  indicated  by  eq.  (12),  the  local  modes  are  normalized  to  carry  the  same  energy  flux  across 
planes  x  =  constant.  We  can  obtain  a  statement  of  energy  conservation  for  a  lossless  range 
dependent  medium  by  substituting  the  local  mode  representation  of  the  stress  displacement 
field  eq.  (13)  into  the  scalar  product  eq.  (11)  and  setting  the  derivative  with  respect  to  x 
equal  to  0.  A  concise  proof  of  this  statement  of  energy  conservation  for  the  general  3-D 
problem  has  been  given  by  Tromp  (1994).  Carrying  out  the  indicated  differentiation  and 
employing  eq.  (14)  yields 


dx 


d 


dx 

E 

9 


'dcgix)  ^  ,  dc*{xy 
-c{x)  +  Cg{x)-^ — 


dx  ^  dx 

E  +  B*^)  Cr{x)cl{x) 


q,r 

0. 


We  have  used  the  fact  that 


(21) 


(E-®«f<w)  =  E'^W  (E-s^Jw)  ■  P2) 

q  \  r  Jr  \  q  / 

since  the  mode  indices  {q,  r}  are  summed  over  the  same  set. 

The  only  way  for  eq.  (21)  to  be  satisfied  generally  is  for  the  coupling  matrix  to  be  anti- 
Hermitian,  i.e. 


=  -b;,.  (23) 

This  anti-Hermiticity  is  a  necessary  consequence  of  energy  conservation  in  a  lossless  medium. 
It  can  be  seen  by  inspection  that  the  coupling  matrix  Bgr  given  by  eq.  (Al)  is  anti-Hermitian. 
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5  THE  COUPLED  ENERGY  EQUATIONS 


As  the  frequency  increases,  the  discrete  mode  spacing  becomes  fine  enough  that  the 
displacement-stress  field  may  be  treated  as  a  continuum  of  trapped  modes.  Random  fluctu¬ 
ations  of  material  parameters  and  layer  thicknesses  and  boundary  slopes  can  cause  energy 
to  be  exchanged  among  modes.  Within  the  limits  of  a  very  well  defined  set  of  assumptions, 
it  is  possible  to  derive  a  diffusion  equation  to  describe  this  energy  exchange.  If  the  medium 
consists  of  densely  packed  scatterers  imbedded  in  a  matrix,  and  if  the  field  is  characterized  by 
a  large  number  of  closely  spaced  modes,  all  the  propagating  energy  is  scattered  energy.  This 
leads  us  to  relinquish  information  about  individual  mode  phases  in  exchange  for  a  descrip¬ 
tion  of  the  average  energy  transport.  Marcuse  (1974)  derived  coupled  power  equations  for 
propagation  in  an  optical  waveguide  with  random  rough  boundaries.  While  our  derivation  of 
the  coupled  energy  equations  roughly  parallels  that  of  Marcuse  (1974),  there  are  differences. 
Marcuse  (1974)  assumed  that  the  coupling  matrix  Bgr  was  independent  of  range,  which  is 
not  necessary.  Another  difference  is  that  Marcuse’s  final  results  for  the  energy  diffusivity 
depend  on  the  autocorrelation  of  the  boundary  function  itself.  Our  results  depend  on  an  au¬ 
tocorrelation  functional  on  the  horizontal  gradients  of  the  boundary  functions  and  material 
parameters. 

Referring  to  eq.  (21),  the  average  energy  flux  of  local  mode  q  across  planes  x  =  constant 
is 


E,  =  {c,<} .  (24) 

The  angle  brackets  (•)  indicate  an  ensemble  average  taken  over  statistically  similar  realiza¬ 
tions  of  the  medium. 

Differentiating  eq.  (24)  with  respect  to  a:,  we  get 
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and  substitute  eq.  (14)  to  get 


del 

^‘^dx 


(26) 


dx 


+ 


=  ^(^Bgr{x)Cr{x)cl{x)'^  +CC,  (26) 

where  cc  refers  to  the  complex  conjugate  of  the  first  expression  on  the  right  hand  side  of 
eq.  (26). 

The  displacement-stress  field  in  a  random  medium  characterized  by  heterogeneities  dis¬ 
tributed  over  a  range  of  scales  is  statistically  nonuniform.  We  deal  with  this  statistical 
nonuniformity  by  dividing  the  medium  into  a  sequence  of  approximately  statistically  uni¬ 
form  segments  of  length 


C{(x>)  —  X  —  x'.  (27) 

for  a  specified  frequency  band  (Wu  &  Aki  1985b).  We  also  assume  that  the  frequency 
dependent  correlation  length  scale  of  the  medium  heterogeneities  l{u))  is  small  compared  to 
the  statistically  uniform  region  (Fig.  1) 

l{uj)  <  C{lo).  (28) 

We  assume  that  over  the  region  of  approximate  statistical  uniformity  JC,{u}),  we  can  ap¬ 
proximate  the  forward  scattered  mode  amplitudes  Cq{x)  and  Cr{x)  with  the  first  two  terms 
from  the  power  series  expansion  of  the  matrix  exponential  eq.  (17).  So  using  the  first  two 
terms  of  eq.  (17)  substitute 


N 

Cr{x)  ^  Cr(x')  +  ^  Arn(x)Cn(x')  (29) 

71=0 

and 
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(30) 


»  <(x')  +  Y, 

r=0 

for  Cr  and  c*  in  the  first  term  of  eq.  (26).  Corresponding  substitutions  are  made  in  the 
complex  conjugate  term,  yielding 

+  Y^(^Bgr{x)Arn{x)Cn{x')c*{x')'j 
r^n 

+  Y.{B<,r{x)AlXx)Cr{x')cl{x')) 

r 

+  Yj  {Bqr{x)Arn{x)A*^{x)Cn{x')cl{x')'^ 

r^n 

+  cc.  (31) 

Because  of  the  double  summation  over  r  and  n,  the  elements  of  the  product  terms  in  eq.  (31) 
commute  without  requiring  the  interchange  of  indices. 

By  making  several  assumptions  about  the  statistics  of  the  random  heterogeneities,  eq.  (31) 
can  be  reduced  in  a  systematic  way  to  yield  the  coupled  energy  equations.  We  assume  that 
the  range  dependent  displacement-stress  field  amplitudes  Ci(x)  are  slowly  varying  over  the 
correlation  distance  l(uj)  of  the  medium  properties.  This  permits  us  to  express  the  ensemble 
averages  of  the  products  in  eq.  (31)  as  products  of  the  ensemble  averages  (Bendat  &:Piersol 
1986,  pp.  465-471) 


+  Y  (Bgr(x)Arn(x)}  (^Cn(x')c*(x')j 

r,n 

+  Y  (Bgr(x)A*^(x))  (Cr(x')c;(x')) 

r 

+  Y  {Bqr{x)Arn{x)Aqr{xy)  {Cn{x')cl{x')) 

r,n 
+  CC. 


(32) 
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The  coupling  matrix  terms  Bqr,  and  the  horizontally  integrated  matrices  Am  and 
are  sums  of  terms  directly  proportional  to  the  horizontal  gradients  of  the  density  [dp/dx), 
stiffness  tensor  (dcijki/dx)  and  layer  boundary  configurations  (dh/dx).  We  assume  that 
the  distribution  of  horizontal  gradients  of  these  medium  properties  within  the  ensemble  of 
statistically  similar  models  are  zero  mean  so  that 

(Bqr)  =  (Am)  =  (Aqr)  =  0.  (33) 

The  zero  mean  assumption  allows  us  to  eliminate  the  first  term  on  the  right  hand  side  of 
eq.  (32).  We  have  already  assumed  that  the  coupling  is  weak  enough  to  approximate  the 
matrix  exponential  eq.  (17)  by  the  first  two  terms  of  its  power  series.  The  third  term  on  the 
right  hand  side  of  eq.  (32)  will  therefore  be  much  smaller  than  the  second  and  third  terms, 
and  we  ignore  it. 

The  second  term  in  eq.  (32)  contains  a  double  summation  over  r  and  n.  Physically,  this 
term  includes  the  effects  of  all  modes  n  which  couple  to  all  modes  r  which  then  couple  to 
mode  q.  We  make  the  assumption  that  the  main  contribution  to  the  double  summation 
comes  from  modes  n  =  q.  This  is  consistent  with  the  previous  assumption  that  the  mode 
amplitudes  can  be  expressed  by  the  first  two  terms  of  the  series  expansion  of  the  matrix 
exponential. 

After  dropping  the  first  and  fourth  terms  on  the  right  hand  side  of  eq.  (32),  eliminating 
the  summations  over  n  in  the  second  and  third  terms,  and  using  the  anti-Hermiticity  of  Bqr, 
we  have 


-  'E{B,rA;r){cq{x')cl{x')) 

+  'E{BgrA*qr){Cr{x')c;{x')) 

r 

+  cc.  (34) 
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The  remaining  terms  proportional  to  (^Bgr{x)A*j.{x)^  can  be  rewritten  using  the  definition 
of  Aqr  (eq.  (20)),  so  that  the  coupling  terms  form  the  autocorrelation  of  the  coupling  matrix 


(b^(x)A;,{x))  =  {B^(x)Bl,(x"))  dx" 

=  f  'Jlqr{r)  dr,  (35) 

J  —  OO 

where 

K„{t)  =  (B^(x)B’^{x  -  t))  (36) 

is  the  autocorrelation  of  the  coupling  matrix,  and 


r  =  rr  —  x" .  (37) 

The  lower  limit  of  the  integral  has  been  extended  from  x'  to  —  oo  because  we  assume  that 
TZqrir)  contributes  only  over  the  distance  of  a  medium  correlation  length  l{u)),  and 


l{u))  X  —  x'.  (38) 

Since  Bqr  is  evaluated  at  x  and  the  integration  for  Agr  is  over  x",  we  are  able  to  take  Bgr 
inside  the  integral,  and  introduce  the  correlation  variable  t.  The  complex  conjugate  term 
contributes  the  complex  conjugate  of  the  integral  term  we  already  have.  The  final  result  for 
the  coupled  energy  equation  is  then 

BE 

-^  =  ^bqr{E,~Eq)  (39) 

where 

/+c» 

Tlgrix)  dr,  (40) 

-OO 
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the  symmetric  energy  coupling  matrix,  is  essentially  the  spatial  autospectrum  of  the  square 
of  the  coupling  matrix  Bgr-  It  can  be  seen  from  eq.  (Al)  that  the  leading  order  phase 
dependence  of  the  integrand  of  eq.  (40)  will  be  exp{i{k^  —  k'^)r)  which  gives  bgr  the  form  of 
an  autospectrum  to  leading  order. 

Referring  to  eq.  (Al),  it  can  be  seen  that  Bgr  is  a  sum  of  terms  made  up  of  field  com¬ 
ponents,  the  vertical  derivatives  of  the  field  components  and  horizontal  derivatives  of  the 
field  components.  If  we  write  Bgr  to  emphasize  the  horizontal  derivatives  of  the  material 
properties. 


where  f{cijki{x))  is  used  to  indicate  any  of  the  Qij  or  matrix  products  of  the  Cij,  and  if 
we  can  assume  that  B^]),  B^gj)  and  Bfj)  are  constant  over  the  statistically  uniform  regions 
C{lo),  we  find  that  bgr  is  proportional  to  a  sum  of  terms  consisting  of  the  spatial  autospec¬ 
tra  and  cross-spectra  of  dp{x)/dx,  df{cijki{x))/dx  and  dh{x)/dx.  If  we  are  able  to  make 
this  approximation,  specification  of  the  random  properties  of  the  medium  is  reduced  to  the 
specification  of  the  corresponding  autospectra  and  cross-spectra  of  the  medium  properties. 
In  addition  the  coupling  matrix  will  only  have  to  be  computed  once  for  every  statistically 
uniform  segment. 


6  THE  ENERGY  DIFFUSION  EQUATION 

Eq.  (39)  is  an  infinite  set  of  coupled  first  order  equations  that  describes  the  evolution  of 
mode  energy  in  a  random  medium  with  a  correlation  function  that  can  be  computed  once 
the  form  of  the  heterogeneities  has  been  specified.  Eq.  (39)  is  not  particularly  useful  in  the 
form  given  because  the  eigenvalues  of  the  local  modes  in  the  random  medium  must  be  found, 
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the  eigenfunctions  must  be  calculated  and  the  range  dependent  energy  coupling  matrix  must 
be  computed  before  eq.  (39)  can  be  solved.  By  employing  several  additional  assumptions, 
we  can  reduce  eq.  (39)  to  a  diffusion  equation  for  the  diffusion  of  energy  in  the  random 
medium.  Our  derivation  parallels  that  of  Gloge  (1972),  who  derived  a  diffusion  equation  for 
the  propagation  of  energy  in  a  dielectric  optical  waveguide  with  random  index  of  refraction 
perturbations. 

The  elements  of  the  coupling  matrix  eq.  (Al)  depend  on  the  inverse  of  the  eigenwavenum- 
ber  spacing.  Tromp  (1994)  has  shown  that  the  the  relative  coupling  strength  between  modal 
overtones  is  inversely  proportional  to  the  square  of  the  local  wavenumber  spacing 

—  (x\kq-  kr\~‘^  (42) 

Cq 

so  the  strength  of  the  energy  coupling  is  therefore 

^  a  1^,  -  kr\-\  (43) 

kjq 

As  the  frequency  increases,  the  eigenwavenumber  spacing  decreases.  The  {Ak)~'^  depen¬ 
dence  of  the  energy  coupling  indicates  that  the  coupling  interaction  between  immediately 
neighboring  modes  will  be  much  stronger  than  more  distant  intermode  interactions.  If  we 
assume  that  nearest  neighbor  coupling  between  modes  comprises  the  dominant  intermode 
interactions,  eq.  (39)  can  be  written  as  the  differential-difference  equation 

^  -  E,)  +  -  E,).  (44) 

As  the  frequency  increases,  the  mode  spacing  becomes  smaller  and  smaller.  In  the  limit 
in  which  the  mode  spacing  goes  to  zero,  the  mode  spectrum  is  represented  by  a  continuum, 
and  the  differences  in  eq.  (44)  can  be  replaced  with  differentials 


Eq+,-Eq^^Aq. 


dq 


(45) 
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This  allows  us  to  rewrite  eq.  (44)  as 


dx 


—  Eq)  bq^q-i{Eq  -E'g-l) 

=  A,k  ^ 


"'9+l>9' 


dq 


bq,q-l' 


dq 


Finally  we  collapse  the  final  difference  in  the  same  way  to  get 


(46) 


dE„ 


a  (  as; 


(47) 


dx  \dq  dq  ) ) 

We  have  dropped  the  second  subscript  on  bqr,  and  set  Aq  =  1. 

The  mode  number  q  is  not  particularly  useful  as  the  continuous  independent  variable. 
However  at  this  point  we  have  quite  a  bit  of  flexibility  for  the  choice  of  continuous  inde¬ 
pendent  variable  to  use  as  we  take  the  continuum  limit.  Possible  choices  include  the  local 
horizontal  wavenumber  kq,  the  local  horizontal  slowness  in  the  propagation  direction  Px,q, 
the  local  phase  velocity  or  the  local  ray  equivalent  angle  6q.  A  final  choice  would  be 
made  according  to  application,  mathematical  perspective  one  would  like  to  take  and  ease  of 
numerical  implementation.  In  any  case  we  need  to  be  able  to  relate  the  mode  number  q  to 
the  continuous  independent  variable  of  choice.  If  we  let  A  represent  any  appropriate  spectral 
variable  we  may  want  to  choose,  then  the  local  dispersion  relation 


nK<l)  =  0  (48) 

provides  an  implicit  relation  between  the  mode  index  q  and  A.  Differentiation  of  eq.  (48) 
provides  the  needed  relation  to  transform  the  derivatives  with  respect  to  q  in  eq.  (47)  to 
derivatives  with  respect  to  A. 

Our  final  result  is  then 
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where 


dx 


Q  \qA,x) 


A,a:)fe(A,2;) 


dE{X,x)\\ 

dX  j/’ 


(49) 


Q{q,X,x) 


dX 


is  the  local  mode  density  function,  and 


(50) 


/+C» 

7?.(r)  dr, 

-OO 


(51) 


which  is  the  continuum  limit  of  eq.  (40)  described  below. 

The  form  of  b{X,x)  follows  from  the  continuum  limit  of  eq.  (36)  and  eq.  (Al).  The 
coupling  matrix  Bgr  (eq.  (Al))  is  a  function  of  the  difference  of  the  eigenwavenumbers  Ak  = 


“  Ak{x) 


Bgr{x)e 


iip{Ak,x) 


(52) 


where  Bqr  is  the  part  of  Bgr  that  does  not  depend  on  Ak  and  ip{Ak,  x)  is  the  phase.  We  can 
express  Ak{x)  in  terms  of  the  mode  density  function  Q{q,  X,  x)  as 


A  ,  ^  dXdk 

.  dk  ,  , 

=  (9,A,a;). 


(53) 


Now  rewrite  the  autocorrelation  function  eq.  (36)  as 


{Bgr{x)B*  {X  -  t))  = 


{Aqy 


;{B{x)B{x  -  r)), 


where 


(54) 
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B{x)  =  oiQ,  A,  x)e{x)e<(^««='»  (55) 

and 

B*{x-t)  =  ^(9,A,a;-r)B*(x-r)e-'(^'?’^("-^».  (56) 

The  tildes  over  '^{x)  and  'ij){x  —  t)  indicate  that  the  A  A:  term  in  the  phases  has  also  been 
replaced  by  the  expression  given  in  eq.  (53).  Allowing  the  mode  spacing  to  approach  zero, 
and  setting  A5  =  1  in  eq.  (54),  we  get 

/+00  /■+00  _ 

Tlqr{T)dT  ->  /  {B{x)B*  {x  -  t)) 

/+00 

'R-{t)  dr 

-00 

=  b{X,x).  (57) 

If  it  occurs  that 

Q{q,  A,  x)  =  Q{q,  X,x  —  t)  ^  constant,  (58) 

for  all  X  and  A,  then  can  be  factored  from  b{X,x)  and  canceled  by  the  inverse  factors 
on  the  right  hand  side  of  eq.  (49)  to  eliminate  the  mode  density  from  the  energy  diffusion 
equation.  Also  note  that  if  we  choose  X  =  k  then  dX/dk  =  1,  and  if  we  choose  X  —  Px  then 
dX/dk  —  l/cj.  Either  of  these  choices  for  A  will  significantly  simplify  the  computation  of 
b(X,  x). 

Eq.  (49)  is  a  parabolic  partial  differential  equation  that  describes  the  range  dependent 
diffusion  of  energy  in  A-space  in  a  perfectly  elastic  2-D  random  medium  that  is  so  strongly 
scattering  that  all  the  energy  is  scattered  energy.  The  diffusivity  b{X,x)  is  a  range  and 
frequency  dependent  continuous  functional  of  the  displacement-stress  field  and  the  range 
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gradients  of  the  medium  properties.  Because  we  have  taken  the  continuum  limit  in  the 
derivation  of  eq.  (49),  the  diffusivity  b(X,x)  can  be  computed  using  any  program,  e.g.  re¬ 
flectivity  or  finite  difference,  that  can  solve  for  the  displacement-stress  field  as  a  function 
of  depth  in  a  plane  layered  medium.  The  field  quantities  in  the  coupling  matrix  Bqr  from 
which  the  diffusivity  6(A,  x)  is  derived,  are  the  field  quantities  for  a  plane  layered  medium. 
The  range  dependence  enters  through  the  range  gradients  of  the  medium. 

In  order  to  solve  eq.  (49),  we  must  be  able  to  evaluate  g{q,\,x).  There  are  several 
special  cases  for  which  Q{q,  A,  x)  can  be  evaluated  analytically  subject  to  some  additional 
approximations.  From  a  practical  numerical  standpoint,  we  need  to  know  the  density  of 
modes  in  the  region  of  interest  in  the  spectral  space  of  choice  {k,  p^,  0,  etc.).  This  is 

straightforward  for  SH  waves  and  acoustic  waves  which  are  governed  by  a  second  order  scalar 
differential  equation  that  defines  a  Sturm-Liouville  problem.  The  eigenvalues  are  predictably 
and  evenly  spaced,  and  it  is  relatively  straightforward  to  determine  the  density  of  eigenvalues 
in  some  region  of  spectral  space.  Determining  the  density  of  eigenvalues  for  the  higher  order 
P-SV  problem  is  more  difficult  because  of  their  irregular  spacing. 

However,  Woodhouse  (1981,  1988)  has  been  able  to  construct  a  function  from  the  minors 
of  the  propagator  matrix  that  exactly  yields  the  density  of  modes  in  a  specified  region  of 
spectral  space  for  the  general  P-SV  problem.  The  evaluation  of  this  function  requires  only 
two  integrations  of  the  system  of  differential  equations  for  the  P-SV  problem.  The  existence 
and  form  of  this  function  were  conjectured  quite  some  time  ago  (Woodhouse  1981),  and 
a  complete  proof  was  given  just  recently  (Woodhouse  1988).  An  effective  algorithm  for 
computing  the  mode  density  function  has  been  given  by  Gomberg  &  Masters  (1988).  We 
discuss  the  more  general  P-SV  case  first,  and  then  the  cases  for  which  it  is  possible  to  find 
an  approximate  analytical  expression  for  0{q,X,x). 

Woodhouse’s  (1988)  function  is 
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(59) 


which  is  the  number  of  eigenvalues  (modes)  of  an  n  x  n  system  of  equations  between  the 
values  Aa  <  Afc.  We  take  as  the  average  mode  density 

Q{q,  A,  x)  ^  y  ■  (60) 

Eq.  (60)  should  provide  a  reasonable  approximation  to  the  average  mode  density  at  high 
frequencies  where  the  modes  are  spaced  closely  together.  It  is  not  necessary  to  precisely  locate 
eigenvalues,  nor  is  it  necessary  to  construct  eigenfunctions  in  order  to  compute  eq.  (60).  It 
is  only  necessary  to  count  the  number  of  zero  crossings  of  a  certain  product  of  minors  of  the 
propagator  matrix.  The  mode  density  Q{q,  A,  a:)  is  also  range  dependent,  and  must  be  pre¬ 
computed  locally  and  stored  along  with  the  quantities  necessary  to  construct  the  diffusivity 
functional  6(A,  x). 

Eq.  (60)  is  quite  general,  and  valid  for  an  arbitrary  fluid-elastic  medium.  We  should 
be  able  to  construct  a  numerical  solution  to  the  energy  diffusion  equation  (49)  for  any  2-D 
fluid-elastic  random  medium  of  interest  that  satisfies  the  assumptions  already  stated.  With 
some  additional  assumptions  it  is  possible  to  obtain  a  precise  analytical  expression  for  the 
mode  density. 

We  choose  the  horizontal  slowness  as  our  independent  spectral  variable,  and  restrict 
it  to  lie  in  the  range  <  Px  <  Pmin  where  and  ^min  are  the  minimum  compressional 
and  shear  speeds,  respectively.  In  this  slowness  range  the  shear  waves  are  propagating  and 
the  compressional  waves  are  evanescent.  If  the  surface  layers  consist  of  water  overlying  very 
low  shear  speed  sediments,  we  take  amin  equal  to  the  minimum  sound  speed  in  the  water, 
and  take  (5min  equal  to  the  shear  speed  in  the  sediments  at  the  water-sediment  interface.  We 
assume  further  that  the  sediment  shear  speed  is  less  than  the  compressional  speed  in  the 
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water.  If  the  medium  can  be  modeled  as  piecewise  smooth,  the  high  frequency  asymptotic 
limit  of  the  dispersion  relation  is  (Kennett  1983,  p.297) 

2ojt^{p^)  ~  {2q  +  1/2)%  +  (61) 

where 

fZp{px)  . 

MPx)^  l4{z)dz,  (62) 

«/  0 

the  intercept  time,  is  the  integral  of  the  vertical  shear  wave  slowness  p^^{z)  from  the  surface 
to  the  shear  wave  turning  depth  Zji(px),  and  'i>(LO,Px)  is  a  frequency  dependent  phase  shift 
from  the  effect  of  reflection  from  all  interfaces,  including  the  free  surface.  If  the  only  interface 
is  a  solid  free  surface  then 


(63) 


where  R^^{px)  is  the  SS  free  surface  reflection  coefficient,  which  is  frequency  independent. 
Prom  eq.  (61)  we  get 


dq 


dpx  %  dpx 

The  first  term  inside  the  braces  on  the  right  hand  side  can  be  simplified  further. 


(64) 


where 


duJTp 

dpx 


=  ^Tp{p^)Us{Px) 


U}Xs{Px), 


(65) 


Vs=^-^ 

udpx 

is  the  local  shear  wave  group  velocity,  and 


(66) 
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(67) 


Xs(p,)  =  Ts(p.)Us(p.),  (69) 

we  have  for  the  mode  density  from  eq.  (64) 

q{q,Px,x)  =  ^  =  -i  {ujp^Xs{Px)Us{jPx)  +  •  ('^0) 

OPx  ^  I  "  ^Px  ) 

By  combining  eq.  (70)  with  eq.  (49),  we  have  an  equation  that  describes  the  range 
dependent  diffusion  of  energy  in  slowness  space  for  the  slowness  range  indicated  above, 

=  ?-■(«, p.,x)  ||-  (71) 

where  Q{q,Px,x)  is  given  by  eq.  (70),  and  from  eq.  (57) 

„  r+oo 

b{Px,  x)  =  /  K{t)  dr.  (72) 

J—OO 

The  factor  of  is  from  (dpx/dkY  in  eq.  (57),  which  we  have  removed  from  TI{t).  All  the 
medium  dependent  factors  are  computed  locally  with  any  suitable  one  dimensional  code  that 
yields  the  complete  displacement-stress  field  as  a  function  of  depth.  The  range  dependence 
enters  explicitly  when  the  diffusivity  functional  b{px,x)  is  constructed. 

Energy  diffusion  of  SH  waves  or  pure  acoustic  waves  is  also  governed  by  equations  of 
the  form  of  eqs.  (70)  and  (71).  The  only  difference  is  that  the  free  surface  contribution  to 
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'i>{LO,Px)  will  be  zero.  This  is  because  there  is  no  wave  type  conversion  for  SH,  or  for  pure 
acoustic  waves,  in  2-D  media.  If  the  medium  is  slowly  varying,  and  the  free  surface  is  the 
only  interface,  then 


^{u>,Px)  =  0. 


(73) 


If  we  expand  the  derivatives  with  respect  to  in  eq.  (71),  we  get  an  equation  the  form 
of  a  convection-diffusion  equation  with  non-constant  coefficients 


dx 


,2/  .d‘^E{px,x)  2/„  ^^dE{px,x) 


E>\px,x)- 


dpl 


r  {Px,x)- 


dpx 


where  the  diffusivity 


(74) 


'D'‘{Px,x)  =  0  '^{q,Pi,x)h(ji„x),  (75) 

and  T‘^{px,x)  characterizes  the  rate  at  which  energy  spreads  from  a  particular  region  of 
Pa;-space 


=  -Q~^{q,Px,x)—  [Q~^{q,Px,x)h{px,x)Y 


dpx 


(76) 


Finally  by  introducing  the  transformation 


E{px,x) 


=  S{px,x){r{px,x)} 


(77) 


the  first  derivative  with  respect  to  Px  can  be  eliminated  from  eq.  (74),  reducing  it  to  a 
Schrodinger  equation 


d'^e 

dpi 


-f  V{px,  x)8  =  2?  ^ 


dx^ 


(78) 
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with  the  potential 


= lA  av  _  1  fiv  _  AA  r  iAA)  Vv 

’  '  2dp^\v)  i\v)  2V^dxJ  \V{if^,x)) 

The  Schrodinger  equation  form  may  be  useful  because  of  the  large  amount  of  work  that  has 
gone  into  its  numerical  solution.  A  drawback  may  be  the  complicated  form  of  the  potential 
V{pxix)^  and  the  boundary  conditions  discussed  in  the  next  section. 


7  INITIAL  AND  BOUNDARY  CONDITIONS 

Eq.  (71)  is  first  order  in  x  and  second  order  in  p^.  In  order  to  completely  define  the 
diffusion  problem,  we  need  an  initial  condition  and  two  boundary  conditions.  As  an  initial 
condition,  which  does  not  have  to  satisfy  the  boundary  conditions,  we  specify  the  source 
spectrum  at  a;  =  0 


£/(p,,0)=<S(p,).  (80) 

For  the  particular  problem  defined  above  at  the  end  of  Section  6,  the  boundary  conditions 
are  to  be  specified  at  the  endpoints  of  a  finite  region  of  spectral  space  <  Px  <  ^min- 
Within  this  region  of  Px  —  space  the  propagating  waves  are  almost  pure  SV.  This  spectral 
region  was  chosen  as  an  illustrative  example  for  which  the  mode  density  function  was  easy  to 
calculate.  First  we  will  specify  the  appropriate  boundary  conditions  for  this  specific  problem, 
then  comment  on  how  they  must  be  changed  for  more  general  problems. 

For  Px  >  Ihe  compressional  waves  are  evanescent  and  carry  no  energy.  For  px  < 
energy  in  the  SV  waves  is  radiated  away  as  P  waves.  We  assume  that  this  energy  is  carried 
away  and  completely  lost  to  the  system.  Because  of  this  energy  loss  we  choose 

^(«mLA)=0  (81) 
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for  the  lower  boundary  condition,  which  corresponds  to  infinite  loss. 

For  the  upper  boundary  condition  at  we  also  choose 

=  (82) 

The  reason  for  this  choice  is  that  for  values  of  the  propagating  SV  waves  can  lose 

energy  to  the  fundamental  Rayleigh  wave  or  low  speed  Stoneley  waves.  This  choice  for  the 
upper  boundary  condition  is  somewhat  artificial  since  energy  can  flow  into  Rayleigh  waves, 
which  are  ignored  by  this  choice. 

If  we  give  up  being  able  to  analytically  calculate  the  mode  density  function  in  the  interests 
of  covering  the  complete  slowness  range  0  <  Pa,  <  oo,  we  need  to  reexamine  the  boundary 
conditions.  For  the  lower  boundary  condition  at  Pmin  —  0,  we  must  choose 


dE 

dpx 


=  0,  for  Pa;  =  0. 


(83) 


Because  we  have  assumed  that  there  is  no  backscattering,  the  energy  flux  into  regions  of 
spectral  space  for  which  Pa;  <  0  must  be  zero.  Eq.  (83)  is  the  mathematical  statement  of  the 
zero  flux  condition. 

We  must  choose  a  finite  upper  limit  for  the  upper  boundary  since  there  are  no  propagating 
waves  for  p^  >  Pmax,  where  Pmax  is  the  slowness  for  the  fundamental  Rayleigh  or  Stoneley 
wave,  whichever  has  the  greatest  slowness.  Therefore,  at  the  upper  boundary  we  also  have 
a  zero  flux  condition 


dE 

dpx 


0,  for  Px  —  Pmaxi 


(84) 


indicating  that  no  energy  is  allowed  to  flow  into  the  spectral  region  for  which  Px  > 


Pmax- 
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If  the  Schrodinger  equation  form  eq.  (78)  is  employed,  the  zero  flux  boundary  conditions 


eq.  (83)  and  eq.  (84)  assume  the  form 


dpa:  ^  \T)) 


0. 


(85) 


Mixed  boundary  conditions  of  the  form  eq.  (85)  complicate  the  solution  process  in  general. 
The  zero  energy  boundary  conditions  such  as  eq.  (81)  and  eq.  (82)  retain  their  same  form, 
as  does  the  initial  condition  eq.  (80). 


8  DISCUSSION  OF  ELEMENTARY  SOLUTIONS  TO  THE  DIFFUSION 
EQUATION 

The  energy  diffusion  problems  described  by  eq.  (49)  or  eq.  (71)  and  the  boundary  conditions 
given  in  the  previous  section  must  be  solved  numerically  for  any  realistic  problem  of  interest. 
However  some  qualitative  insight  can  be  gained  by  examining  simple  solutions  for  constant 
diffusivity.  We  define  two  example  problems  that  illustrate  the  points.  Because  we  specify 
constant  diffusivity,  the  “convection”  term  in  eq.  (74)  vanishes  since  the  coefficient  of  that 
term  is  the  derivative  of  the  diffusivity. 

Problem  I. 

The  first  example  problem  is 


dE 

—  =  TE — - 
dx  dpi 


X>Q  and  Pmin  <P<  Pmax, 


(86) 


with  initial  condition  eq.  (80)  and  boundary  conditions 


E{pmint  x) 


dE{pirriaxj  X^  _  q 

dpx 


(87) 
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The  constant  diffusivity  is  denoted  by  V^,  and  has  units  of  [slowness‘s /length].  Problem  I 
assumes  zero  energy  flux  above  a  maximum  value  of  =  Pmax  and  assumes  that  any  energy 
scattered  to  compressional  waves  at  p^  —  Pmin  is  completely  lost  to  the  system.  This  problem 
defines  an  eigenvalue  problem  for  the  energy  E{px,  x)  in  the  region  <  P  <  Pmax  that 
can  be  solved  by  a  number  of  different  methods,  It  is  treated  in  detail  by  most  textbooks 
on  partial  differential  equations,  so  we  simply  state  the  Fourier  series  representation  of  the 
solution 


OO 

E{px,  x)  =  Y,  •  {tan(A„p^a^)  sin(A„p^)  +  cos(A„p^)}  ,  (88) 

n—l 

where  the  eigenvalue  A„,  which  is  not  related  in  any  way  to  the  spectral  parameter  A  used 
earlier,  is 


Aji  — 


7r(2n  —  1) 


2(Pmax  Pmin) 

and  n  is  the  summation  index.  The  expansion  coefficients  are 


(89) 


{(» 


—  2  ^  {Pmax  Pmin)  ^  +  Pmin 


^  rPmax 

/ 

■'  ''Pmin 


S{p^)  cos{XnPx)  dp^.  (90) 


There  are  several  things  to  notice  about  the  form  of  the  solution  eq.  (88) .  It  is  apparent 
that  as  the  propagation  distance  in  the  strongly  scattering  medium  x  oo  that  the  energy 
E  ^  0.  All  the  energy  is  eventually  scattered  away  to  body  waves  and  lost  to  the  system. 
Also  any  structure  that  may  initially  exist  in  the  Pa;-spectrum  is  smeared  out  as  the  prop¬ 
agation  distance  increases.  The  sharpest  peaks  in  the  spectrum  correspond  to  the  largest 
eigenvalues  of  the  solution  eq.  (88),  but  the  spectral  components  with  the  largest  eigenvalues 
suffer  the  largest  damping  proportional  to  exp(— n^x)  with  range.  The  effect  of  a  strongly 
scattering  medium  is  to  homogenize  the  spectrum  and  destroy  any  coherent  signal. 
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Problem  II. 


The  second  example  problem  is  defined  by  replacing  the  zero  energy  boundary  condition 
at  Px  =  Pmin  by  the  zero  flux  boundary  condition  eq.  (83)  Px  —  0.  Because  of  the  zero 
flux  boundary  conditions  at  both  ends  of  the  spectral  interval,  there  is  no  energy  loss  to  the 
system.  The  solution  is  a  sum  of  a  steady  state  (ss)  part  and  a  transient  (t)  part. 


E{Pxt  —  ^ss  T  £■(, 

where  the  steady  state  part  is  the  total  energy  divided  by  the  interval  length 


(91) 


2  rPmax 

Ess  = - /  S{px)dpx. 

Pmax  -^0 

The  transient  part  of  the  solution  is 


oo 

Et{Px,x)  =  Y.  COs{XnPx), 

n=l 

where 


(92) 


(93) 


and 


HTT 


Pmax 


(94) 


2  rPmax 

an  = -  /  S{px)  cos{XnPx)  dpx-  (95) 

Pmax  ■'0 

In  the  limit  as  x  ->  oo  the  transient  part  of  the  solution  vanishes  and  the  energy  is  distributed 
uniformly  throughout  the  spectral  region  0  <  Pa,  <  Pmax-  As  with  Problem  I  all  coherent 
structure  disappears  from  the  spectrum.  This  is  illustrated  in  Figure  2.  The  remaining 
steady  state  solution  corresponds  to  the  zero  eigenvalue.  For  a  realistic  medium  with  non¬ 
zero  intrinsic  attenuation  there  will  of  course  be  no  steady  state  solution.  The  energy  will 
approach  zero  as  a:  -4  oo. 
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The  only  requirement  for  the  validity  of  the  Fourier  series  representation  of  the  solutions 
to  Problem  /  and  Problem  //is  that  S{px)  be  continuous.  We  have  previously  assumed  in  the 
derivation  of  the  diffusion  equation  that  the  spectrum  could  be  represented  by  a  continuum. 
The  disadvantage  of  the  representations  eq.  (88)  and  eq.  (93)  from  a  computational  stand¬ 
point  is  that  the  series  have  rather  poor  convergence  properties  for  small  x.  For  very  small 
X  the  Green’s  function  representation  of  the  solutions  is  more  efficient.  We  do  not  reproduce 
the  Green’s  functions  for  the  problems  here.  Kevorkian  (1990,  Chapter  1)  has  an  excellent 
discussion  of  the  various  solution  representations  for  the  diffusion  equation  including  a  proof 
of  the  equivalence  of  the  Fourier  series  and  the  Green’s  function  solutions.  A  final  comment 
on  the  Green’s  function  representation  is  that  it  is  very  inefficient  for  large  x  because  a  large 
number  of  image  sources  are  required  in  order  to  satisfy  the  zero  flux  boundary  conditions 
for  either  Problem  /  or  Problem  II. 


9  ATTENUATION 

The  theory  presented  above  is  also  valid  for  arbitrary  dissipative  media  if  the  scalar  product 
eq.  (11)  is  replaced  with  the  more  general  scalar  product  defined  by  Maupin  (1992).  The  form 
of  the  energy  diffusivity  functional  is  such  that  it  can  be  locally  calculated  with  a  reflectivity 
program.  Most  reflectivity  codes  incorporate  arbitrary  attenuation  by  introducing  complex 
frequency  dependent  velocities  or  elastic  moduli.  The  viscoelastic  correspondence  principle 
states  that  the  governing  equations  and  boundary  conditions  for  the  displacement-stress 
fields  are  formally  identical  in  the  perfectly  elastic  and  dissipative  cases.  The  only  difference 
is  that  the  velocities  or  the  elastic  moduli  for  the  perfectly  elastic  media  are  replaced  with 
complex  velocities  or  elastic  moduli  for  the  dissipative  media.  Attenuation  then  enters  the 
problem  in  a  natural  and  straightforward  way,  rather  than  in  the  ad  hoc  manner  used  by 
Gloge  (1972),  Marcuse  (1974)  and  all  the  energy  diffusion  equations  (so  far  as  we  are  aware) 
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originating  from  radiative  transfer  theory. 

One  problematic  practical  issue  is  the  computation  of  the  mode  density  function  A,  a;), 
which  would  have  to  be  done  in  complex  {q,  A)  space.  Woodhouse’s  (1988)  mode  count 
function  u{Xa,  Xb)  for  the  Rayleigh  wave  case  was  derived  for  real  systems  of  equations,  not 
complex  systems.  A  possible  solution  might  be  to  estimate  it  with  its  real  value. 


10  CONCLUSIONS 

We  have  derived  a  very  general  energy  diffusion  equation  (49)  directly  from  the  local  coupled 
mode  representation  of  the  displacement-stress  field  for  a  randomly  heterogeneous  2-D  fluid- 
elastic  medium.  Eq.  (49)  describes  the  diffusion  in  spectral  space  as  a  function  of  range 
in  media  whose  propagation  processes  are  dominated  by  strong  scattering  in  the  forward 
direction.  The  derivation  proceeds  in  a  systematic  way  following  a  well  defined  sequence 
of  assumptions,  and  is  explicitly  frequency  dependent.  Because  of  the  large  number  of 
assumptions  required  to  derive  the  energy  diffusion  equation,  we  summarize  them  here: 

1.  Backscattering  is  neglected.  This  implies  that  our  results  are  most  applicable  to  media 
characterized  by  velocity  heterogeneities  or  very  weak  impedance  heterogeneities. 

2.  Mode  coupling  is  weak  enough  to  express  the  forward  scattered  amplitudes  by  the  first 
two  terms  of  the  series  expansion  for  the  matrix  exponential.  This  approximation  takes 
first  order  mode  coupling  into  account. 

3.  The  statistically  non-uniform  random  medium  can  be  characterized  by  a  series  of  re¬ 
gions  of  frequency  dependent  length  C{u))  that  are  approximately  statistically  uniform. 

4.  The  frequency  dependent  correlation  length  scale  l{uj)  of  the  medium  heterogeneities  is 
much  smaller  than  the  length  scale  of  the  statistically  uniform  regions  (/(cj)  ^  >C(a;)). 
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5.  The  range  dependent  displacement-stress  field  amplitudes  Cq{x)  vary  slowly  over  the 
correlation  distance  l{u!)  of  the  medium. 

6.  The  distributions  of  the  horizontal  gradients  of  the  medium  properties  are  zero  mean. 

7.  The  autocorrelation  of  the  coupling  matrix  'R.qrij)  —  {Bqr(x)B*^{x  —  r))  makes  a 
significant  contribution  only  over  the  distance  of  a  medium  correlation  length  /(w). 

8.  The  mode  spectrum  can  be  represented  by  a  continuum  in  the  high  frequency  limit. 

A  further  useful  approximation  is  that  the  part  of  the  coupling  matrix  Bq,.  that  depends 
directly  on  the  displacement-stress  field  components  is  constant  throughout  the  statistically 
uniform  regions.  If  this  is  the  case,  the  energy  diffusivity  functional  is  seen  to  be  directly 
proportional  to  the  autospectra  and  cross-spectra  of  the  horizontal  derivatives  of  the  medium 
properties. 

The  specific  example  given  for  the  energy  diflFusion  of  Rayleigh  waves  in  the  slowness 
range  Q;“i„  <  Px  <  Pmin  where  and  are  the  minimum  compressional  and  shear 
speeds  of  the  medium,  is  also  formally  identical  to  the  energy  diffusion  equation  for  Love 
waves  or  for  scalar  acoustic  waves  propagating  in  a  waveguide. 

Elementary  solutions  to  the  diffusion  equation  for  relevant  boundary  conditions  clearly 
illustrate  the  the  homogenization  of  the  spectrum  as  the  propagation  distance  increases.  The 
most  singular  features  of  the  initial  spectrum  are  the  most  heavily  damped  with  range. 
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APPENDIX  A:  MAUPIN’S  COUPLING  MATRIX  FOR  FLUID-ELASTIC 
MEDIA 


In  this  appendix  we  include  for  reference,  Maupin’s  (1988)  mode  coupling  matrix  for  a  general 
range  dependent  model  comprised  of  a  fluid  layer  of  depth  hi{x)  lying  over  a  semi-infinite 
solid  halfspace.  The  incompressibility  of  the  fluid  layer  is  designated  by  k,  which  is  a  minor 
deviation  from  Maupin’s  (1988)  notation.  The  designations  hf  and  h'^  indicate  that  those 
terms  are  to  be  evaluated  in  the  fluid  just  above  the  fluid-solid  interface  and  in  the  solid  just 
below  the  fluid-solid  interface,  respectively.  The  dot  (e.g.  p)  indicates  differentiation  with 
respect  to  x.  The  elements  of  the  and  the  Q^j  for  isotropic  and  transversely  isotropic 
media  are  given  in  Maupin  (1988)  and  (1992),  respectively. 
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Figure  Captions 

Figure  1.  Geometry  and  notation  used  for  the  energy  diffusion  problem.  The  randomly 
heterogeneous  region  is  xl  ^  Xr,  is  divided  up  into  N  approximately  statistically  stationary 
regions  Ci{oj),  {i  —  1,N).  Within  any  approximately  stationary  region  jC,i{uj),  we  assume 
where  is  the  frequency  dependent  correlation  length. 

Figure  2.  Schematic  representation  of  the  evolution  of  the  energy  spectrum  with  propaga¬ 
tion  distance.  The  initial  condition  at  x  =  0  must  be  continuous,  but  may  have  sharp  closely 
space  peaks  (a).  The  sharpest  spectral  features  decay  most  rapidly  with  range  according  to 
eqs.  (88)  and  (93)  (b).  In  the  limit  that  x  oo,  the  energy  is  distributed  uniformly  within 
the  propagating  wave  spectral  band  with  energy  amplitude  given  by  eq.  (92)  (c). 
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Introduction 

An  acoustic  signal  propagating  in  shallow  water  can  be  conveniently  represented  as  a 
linear  superposition  of  orthogonal  modes.  Our  definition  of  shallow  water  is  frequency 
dependent  since,  by  "shallow"  we  mean  that  the  water  is  a  "few"  acoustic  wavelengths  in 
depth.  A  "few"  is  somewhat  arbitrary,  but  may  be  as  many  as  100.  The  main  point  is  that 
acoustic  signal  propagation  in  shallow  water  exhibits  interference  and  diffraction 
phenomena  that  are  not  well  modelled  with  ray  theory. 

When  modal  representations  of  the  acoustic  field  have  been  applied  to  shallow  water 
sound  propagation  by  the  ocean  acoustics  community,  it  has  been  fairly  common  practice 
to  include  the  effects  of  attenuation  as  a  first  order  perturbation  to  modal  eigenvalues 
(e.g.  Ingenito,  1973;  Zhou,  1985).  First  order  perturbation  theory  ignores  anelastic  cou¬ 
pling  between  modes  and  requires  that  k/{Q^k)-^\  where  k  is  the  unperturbed 
wavenumber,  hk  is  the  unperturbed  wavenumber  spacing  and  Q  is  the  spatial  quality  fac¬ 
tor.  Because  at  low  frequencies  a  significant  fraction  of  a  shallow  water  acoustic  signal 
propagates  in  low  Q  bottom  sediments,  k/(Qdk)  can  be  O  (1).  This  makes  the  use  of  per¬ 
turbation  theory  invalid,  and  can  introduce  serious  errors  in  mode  sum  acoustic  signal 
synthesis.  Ewing  et  al.  (1992)  report  shear  Q  values  in  the  range  20  -  50  for  continental 
shelf  sediments  off  the  New  Jersey. 

The  severity  of  the  error  resulting  from  the  improper  treatment  of  attenuation  in  mode 
sum  signal  synthesis  has  been  graphically  illustrated  by  Day  et  al.  (1989).  Day  et  al.  cal¬ 
culated  synthetic  seismograms  for  stratified  models  consisting  of  a  high  Q  layer  over  a 
layered  half  space.  The  shear  Q  of  the  underlying  half  space  was  lower  than  the  Q  in  the 
overlying  layer,  and  half  space  shear  speeds  were  lower  than  the  compressional  wave 
speed  in  the  overlying  layer.  The  signals  calculated  from  modal  summation  dramatically 
overestimate  the  effect  of  the  low  shear  Q  on  the  complete  signal.  The  mode  summation 
results  are  compared  with  the  results  from  a  wavenumber  integration  routine  (Apsel,  R.J. 
and  J.E.  Luco,  1983)  that  is  similar  to  SAFARI  (Schmidt  and  Tango,  1986),  at  least  in  the 
way  in  which  attenuation  is  incorporated  into  the  model.  Serious  errors  in  the  mode  sura 
seismograms  are  traced  by  Day  et  al.  to  the  way  in  which  perturbation  theory  is  applied 
to  treat  the  anelastic  problem.  Specifically,  the  difficulty  occurs  with  the  use  of  the 
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unperturbed  elastic  eigenfunctions  in  the  anelastic  problem.  Because  the  problem  is  with 
the  eigenfunctions  themselves,  it  cannot  be  repaired  with  higher  order  perturbation 
corrections  to  the  eigenvalues.  The  problem  does  not  appear  when  using  wavenumber 
integration  routines,  because  the  attenuation  is  incorporated  directly  through  the  use  of 
complex,  frequency  dependent  compressional  and  shear  speeds. 

If  one  wishes  to  retain  the  physical  insight  inherent  in  a  modal  representation  of  the 
acoustic  field,  and  properly  incorporate  the  effects  of  anelasticity,  one  approach  is  to 
invoke  the  correspondence  principle  (e.g.  Leitman  and  Fisher,  1973),  and  solve  for  the 
anelastic  modes  directly.  The  correspondence  principle  states  that  the  equations  of 
motion  for  a  linear  viscoelastic  material  are  just  the  equations  for  a  perfectly  elastic 
material  with  the  elastic  moduli  replaced  with  the  complex,  frequency  dependent  anelas¬ 
tic  moduli.  The  anelastic  moduli  must  be  frequency  dependent  to  preserve  causality. 

The  correspondence  principle  approach  has  been  used  by  Yuen  and  Peltier  (1982)  and 
Buland  et  al.  (1985)  to  model  aspects  of  the  free  oscillations  of  the  whole  earth.  To  a 
limited  extent,  it  has  also  been  applied  to  shallow  water  propagation  problems.  Bucker 
and  Morris  (1965)  employed  the  correspondence  principle  to  solve  for  the  anelastic 
eigenwavenumbers  and  model  the  propagation  loss  for  a  shallow  water  problem  with  a 
fairly  simple  structure. 

We  adopt  a  different  approach  in  that  we  represent  the  anelastic  modes  as  a  complex 
superposition  of  elastic  eigenfunctions.  The  effects  due  to  anelastic  mode  coupling  are 
explicitly  included  and  there  is  no  restriction  on  the  magnitude  of  the  damping.  Our 
approach  is  a  traveling  wave  adaptation  of  Tromp  and  Dahlen  (1990),  who  derived  an 
elegant  solution  for  the  free  oscillations  of  an  anelastic  spherical  earth  in  terms  of  the 
elastic  eigenfunctions  and  eigenfrequencies.  Advantages  of  using  the  elastic  eigenfunc¬ 
tions  as  a  basis  for  the  anelastic  eigenfunctions  are:  1.  The  effect  of  anelasticity  on  indi¬ 
vidual  modes  can  be  examined  in  detail;  2.  The  effect  of  range  dependent  attenuation  can 
be  studied  by  making  the  complex  expansion  coefficients  range  dependent.  If  the 
environment  is  not  geometrically  range  dependent,  we  can  employ  the  same  elastic  basis; 
3.  Used  in  conjunction  with  a  general  range  dependent  coupled  mode  program,  the  pro¬ 
pagation  physics  of  a  strongly  range  dependent  shallow  water  environment  can  be  stu¬ 
died  in  detail.  We  have  the  ability  to  isolate  the  influence  of  the  geometry,  and  different 
aspects  of  the  rheology  of  the  medium  on  a  propagating  shallow  water  signal. 

We  have  derived  a  generalized  eigenvalue  equation  for  the  complex  eigenwavenumbers 
and  complex  coefficients  used  in  the  superposition  of  the  elastic  eigenfunctions  to  con¬ 
struct  the  anelastic  eigenfunctions.  Our  generalized  eigenvalue  equation  is  strictly  linear 
for  the  complex  anelastic  eigenwavenumbers.  This  is  in  contrast  to  the  nonlinear  eigen¬ 
value  equation  for  the  anelastic  eigenfrequencies  of  the  free  oscillations  of  the  earth 
(Dahlen,  1981;  Tromp  and  Dahlen,  1990).  The  reason  for  this  difference  is  our  choice  of 
the  frequency  (O  as  the  independent  variable.  Because  of  the  standing  wave  nature  of  the 
earth  free  oscillation  problem,  co  is  taken  as  the  dependent  variable  in  the  dispersion  rela¬ 
tion.  Since  the  anelastic  moduli  are  frequency  dependent,  the  eigenvalue  equation  for  the 
anelastic  free  oscillations  is  nonlinear. 

Our  derivation  also  includes  the  effects  of  transverse  isotropy,  which  has  a  single  vertical 
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axis  of  symmetry.  It  is  commonly  observed  that  horizontally  propagating  SH  waves  have 
higher  speeds  than  horizontally  propagating,  vertically  polarized  (SV)  or  vertically  pro¬ 
pagating  shear  waves.  In  addition,  higher  compressional  wave  speeds  for  horizontally 
propagating  compressional  waves  than  for  vertically  propagating  compressional  waves 
are  observed  (Ewing  et  al.,  1992).  Berge  et  al.  (1991)  report  shear  wave  anisotropy  of 
12%  - 15%  in  the  same  shelf  region  off  of  New  Jersey  studied  by  Ewing  et  al. 


Definitions  and  notation 

The  derivation  of  the  generalized  eigenvalue  equation  for  the  anelastic 
eigenwavenumbers  will  be  sketched  briefly.  A  more  detailed  treatment  is  currently  being 
prepared  for  publication.  We  adopt  the  notational  conventions  of  Takeuchi  and  Saito 
(1972),  who  treat  seismic  surface  waves  and  free  oscillations  explicitly  for  a  transversely 
isotropic  earth.  As  mentioned  above,  the  case  of  transverse  isotropy  is  probably  the  most 
relevant  for  bottom  interacting  shallow  water  propagation.  Berge  et  al.(1991)  felt  that 
additional  azimuthal  anisotropy  induced  by  ripples  in  the  sediment  surface  or  cracks 
would  be  very  weak.  A  particular  feature  of  transversly  isotropic  media  is  that  the  P  -  SV 
motion  still  decouples  from  the  SH  motion.  This  is  not  true  for  more  general  anisotropy. 
We  include  the  SH  problem  for  completeness,  but  discuss  it  only  briefly.  It  is  the  P  -  SV 
problem  that  describes,  when  the  appropriate  limits  are  taken,  the  propagation  of  an 
acoustic  signal  propagating  in  shallow  water  over  a  plane  layered  viscoelastic  bottom. 

Our  coordinate  system  is  a  right  handed  coordinate  system  with  the  propagation  direction 
along  the  x  axis,  y  positive  into  the  paper,  and  z  positive  upward.  As  mentioned  the  P  - 
SV  (Rayleigh)  motion  decouples  from  the  SH  (Love)  motion.  The  perfectly  elastic  dis¬ 
placements  Ui  and  stresses  O/y  for  P  -  SV  are 

Ux  = 

Uy  =0  (1) 

and 


Ctrr  = 


dy  I 

F-~-kAy3 

dz 


,i{(Ot-kx) 


Oyy  = 


dyi 

F-f^-k(A-2N)y3 

dz 


(2) 

^yz  ~  ^xy  ~  0 
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The  boundary  conditions  on  the  y,-  for  P  -  SV  are 
yi(i  =  1,  2,  3, 4 )  continuous 

^2  =  y4  =  0  at  the  free  surface  z  =  0  (3) 

yi(i  =  1,  2,  3, 4)  ->  0  as  z  ^  -oo 


The  displacements  and  stress  for  SH  motion  are 


and 


Ux  =  u^  =  0 


^yz  ~ 


dyi 

dz 


,  (■  (  tt)f  -  fac ) 


~  ~  ^zz  ~  ^zx  ~  ® 


In  addition,  we  introduce  the  definition  for  y2 

dyi 


y2(z;(0,k)  =  L- 


dz 


so  that 


<^yz  =y2e 


i  { at -kx) 


The  boundary  conditions  for  the  yi  for  SH  are 

yi,y2  continuous 

y2  =  0  at  the  free  surface  z  =  0 

yi,y2  0  as  z 


(4) 

(5) 

(6) 

(7) 

(8) 


In  order  to  minimize  the  notational  overhead,  we  make  no  distinction  between  the  y,  ’s  for 
the  P  -  SV  and  SH  problems.  We  will  concentrate  mainly  on  the  P  -  SV  problem.  The 
meaning  of  the  y/s  will  be  clear  from  context.  The  P  -  SV  and  SH  problems  are  treated 
separately  since  they  do  not  couple  in  transversly  isotropic  media. 

The  above  equations  of  motion  for  a  perfectly  elastic  medium  may  be  represented  in  first 
order  form  as 
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3zbR,L  =  Mr^l  bR,L 


(9) 


where  the  subscripts  R,  L  indicate  whether  we  are  referring  to  the  P  -  SV  [Eq.  (1)  -  (3)] 
or  the  SH  [Eq.  (4)  -  (8)]  systems  of  equations.  The  vectors  bR  l  and  the  matrices  Mr  l 
are  defined  for  P  -  SV  motion  as 

bR  =  (yi,y2,}'3.y4)^  (10) 

and 


Mr  = 


0 


03^  p 


0 


C"* 

0 

0 

-kV/C 


jfcF/C 

0 

0 

a-fVc 


0 

k 


-co^  p 


(11) 


and  for  SH  motion  as 

bL  =  (>'i>3'2)^  (12) 

and 


Ml  = 


0 


it^N-OD^p 


0 


(13) 


In  the  matrices  Mr  r,  w  is  the  real  angular  frequency,  p  is  the  density,  k  is  the  horizontal 
wavenumber  for  a  perfectly  elastic  medium,  and  A,  C,  F,  L  and  N  are  the  five  elastic 
moduli  necessary  to  characterize  a  transversly  isotropic  medium.  When 


A  =  C  =  A,  +  2p  L  =  N  =  |J.  F  =  A, 


(14) 


the  medium  is  isotropic. 

From  this  point  on,  we  drop  the  subscripts  R,  L  on  the  vector  b  and  the  matrix  M.  It  will 
be  apparent  from  context  which  system  we  mean.  The  following  development  will  be  for 
the  P  -  SV  system.  Analogous  results  for  the  SH  system  are  surmnarized  at  the  endof  the 
paper. 

There  are  inherent  symmetries  in  the  equations  of  motion  (Kennett  et  al.,  1978  and 
Thomson  et  al.,  1986)  that  can  be  exploited  to  construct  compact  expressions  useful  for 
very  efficiently  deriving  the  elastic  wave  dispersion  relation,  orthogonality  relationships 
and  other  quantities.  Define  the  matrices  R,  S,  A  and  E  as 
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■  A 

0 

[11 

0 

R  = 

and  S  = 

0 

A 

0  S 

■  0  1  ■ 

■  0 

1  ■ 

A  = 

-1  0 

and 

w  — 

0 

0 

(15) 


(16) 


Using  the  matrices  defined  above,  we  can  form  various  compositions  of  two  stress  fields. 
For  example  for  P  -  SV,  we  can  form 


a,(bTSb)  =  b"^ 


m'^s  +  sm 


b 


(17) 


whereby  we  operate  from  the  left  and  the  right  with  the  same  stress  displacement  vector 
b.  When  both  sides  are  integrated  with  respect  to  z  from  -oo  to  0,  we  obtain  the  disper¬ 
sion  relation  for  Rayleigh  waves  in  a  transversly  isotropic  medium. 


(O^Ij  —  +  ^13  +  I4 


where 


V 

ii  =  j  PCvi^  ^yi^)dz 


I2  =  J  (Lyi^  +  Ay 3^) dz 


yj 

13  =  2/ 


,  dys  _  dyi 

dz  dz 


dz 


\j 

u=  J 


dyi  ^  dy^  ^ 
dz 


dz 


(18) 

(19) 

(20) 

(21) 

(22) 


The  left  hand  side  of  Eq.  (17)  is  a  perfect  differential  in  z,  and  after  integration  it  vanishes 
when  the  boundary  conditions  are  applied. 


Derivation  of  the  generalized  eigenvalue  equation 

The  complex  generalized  eigenvalue  equation  for  the  complex  eigenwavenumbers  is 
derived  in  a  straightforward  manner.  Invoking  the  correspondence  principle,  we 
represent  the  equations  of  motion  for  an  anelastic  transversly  isotropic  medium  as 

=  Me 


(23) 
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where  c  and  M  represent  the  stress-displacement  vector  and  wave  operator  matrix, 
respectively  for  anelastic  media.  For  our  definition  of  M,  we  take 


M  = 


0 

p 


■K 


0 

0 

-kF/C 


K" 


kF/C 

0 

0 

a-fVc 


p 


0 

K 

“-1 

0 


(24) 


The  symbols  A,  C,  F,  L  and  N  are  the  five  complex  frequency  dependent  elastic  moduli 
for  an  anelastic  transversly  isotropic  solid;  K  is  the  eigenwavenumber  for  the  anelastic 
solid;  a  is  the  frequency,  which  we  take  to  be  real  for  propagating  waves. 

We  take 


c  =  c„  (25) 

where  c„  is  an  eigenfunction  of  the  anelastic  medium  and  is  a  solution  of  Eq.  (23).  In 
addition,  we  represent  c„  as  a  complex  linear  superposition  of  the  eigenfunctions  of  the 
perfectly  elastic  problem 


c„ 


^Qm  n^m 

m 


(26) 


is  the  matrix  of  complex  coefficients  that  transforms  the  elastic  eigenfunctions  to 
the  anelastic  solutions. 


Employing  the  matrix  R  defined  above,  we  form  the  composition  of  an  anelastic  mode  c„ 
as  represented  by  Eq.  (26)  and  an  elastic  mode  b„ 


a^Cb^'^Rc^)  =  b/ 


M^R  +  RM 


(27) 


and  integrate  over  z  from  -  oo  to  0.  The  elastic  and  anelastic  problems  satisfy  the  same 
boundary  conditions,  so  the  left  hand  side  of  Eq.  (27)  vanishes  after  the  integration. 

By  assuming  that  the  elastic  eigenfunction  b„  and  the  anelastic  eigenfunction  c„  have  the 
same  real  frequency  so  that 

0)2  =  a2  (28) 

we  arrive  at  the  following  infinite  generalized  quadratic  eigenvalue  equation  for  the  ane¬ 
lastic  eigenwavenumbers  K„ 

Aq„  +  K„Bq„  +  K„2Cq„  =  0 


(29) 
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where 


0 

f.  y 

A  =  - J 

k  2 

A-F^C 

— oo 

>3<”V3 


(m) 


+  k„ 


c 


yi 


(n)y^(m)  _ 


dz 


(30) 


B=  j 


JF 

C 


*(31) 


0  . 


C=J 


A-FVC 


y3^"V3^'”^ 


dz 


(32) 


The  eigenvectors  q„  are  the  columns  of  the  matrix  Q 
Q  =  (...,q„,...) 

By  making  the  assignment 

Ir„  =K„Iq„ 


(33) 


(34) 


the  quadratic  generalized  eigenvalue  problem  Eq.  (29)  can  be  converted  to  a  linear  gen¬ 
eralized  eigenvalue  problem  (Garbow  et  al.,  1977)  at  the  expense  of  doubling  the  dimen¬ 
sion  of  the  system 


■ 

• 

■ 

■ 

A 

B 

qn 

0 

-c 

q« 

0 

I 

r„ 

=  K„ 

1 

0 

r« 

. 

. 

. 

. 

(35) 


Equation  (35)  is  the  main  result  of  this  paper.  The  solution  of  this  linear  generalized 
matrix  eigenvalue  problem  yields  the  complex  eigen3vavenumbers  K„  for  the  modes  of  an 
anelastic  transversly  isotropic  medium  and  the  eigenvectors  q„.  The  eigenvectors  q„  of 
Eq.  (35)  are  the  columns  of  the  transformation  matrix  Qn  m  used  to  construct  the  anelastic 
eigenfunctions  from  the  elastic  eigenfunctions  from  Eq.  (26).  The  linearity  of  Eq.  (35)  is 
an  important  point  and  should  be  contrasted  with  the  result  derived  by  Tromp  and  Dahlen 
(1990)  for  the  free  oscillations  of  the  earth.  Their  equation  (3.5)  is 


Q^+\(Gk) 


(36) 
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The  matrix  Q  is  a  diagonal  matrix  of  eigenfrequencies  for  the  perfectly  elastic  earth; 
is  the  column  of  the  transformation  matrix  Q;  is  the  complex  eigenfrequency  of 
the  anelastic  mode;  and  V(o)  is  an  anelastic  potential  energy  matrix  that  is  a  func¬ 
tional  of  products  of  the  elastic  eigenfunctions  and  the  complex  frequency  dependent 
elastic  moduli.  The  standing  wave  nature  of  the  earth  free  oscillation  problem  leads  to 
the  choice  of  the  frequency  as  the  dependent  variable.  Because  of  the  dependence  of  the 
elastic  moduli  on  frequency,  the  problem,  Eq.  (36),  is  nonlinear.  The  choice  of  frequency 
as  the  independent  variable  for  the  traveling  wave  problem  leads  to  the  linear  problem 
we  have  derived  above,  Eq.(35). 

We  have  also  derived  similar  results  for  the  SH  problem.  We  form  the  compostion 


m'^e  +  hm 


(37) 


As  definitions  of  M  and  b,  we  take  Eqs.  (12)  and  (13)_for  SH  motion.  For  M  we  take  Eq. 
(13)  with  k  replaced  by  K,  and  N  and  L  replaced  by  N  and  L,  respectively.  Likewise  the 
definition  of  c  follows  from  Eq.  (12)  and  (25).  Upon  integrating  Eq.  (37)  with  respect  to  z 
from  0  to  -  oo,  carrying  out  some  additional  algebra  and  again  setting 

0)2  =  ©2  (38) 


we  arrive  at 


k„2l-K„2N-L 


where 


N  = 


-1 


Kf 


and 


L  = 


-1 


dz  dz 


dz 


(39) 


(40) 


(41) 


The  two  terms  6L  and  5N  are  the  complex  frequency  dependent  parts  of  the  two  shear 
moduli  L  and  N.  The  anelasticity  tensor  Cij^i  for  a  linear  viscoelastic  material  can  be 
written 


^ijkl  ~  ^ijkl  ikd)ijkl  (42) 

We  were  able  to  separate  the  perfectly  elastic  part  from  the  frequency  dependent  anelas¬ 
tic  part  for  the  SH  problem.  There  is  no  restriction  on  the  magnitude  of  6L  and  5N.  Also 
note  that  Eq.  (39)  is  linear  in  k„2. 
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A  final  point  is  that  both  Eq.  (35)  and  Eq.  (39)  are  infinite  in  dimension.  Any  practical 
implementation  will  necessarily  employ  a  truncated  mode  set.  This  should  not  be  much 
of  a  restriction,  as  most  low  frequency  shallow  water  signals  can  be  adquately  modelled 
with  a  limited  number  of  modes. 


Conclusions 

The  linear  complex  generalized  eigenvalue  equations  (35)  and  (39)  appear  to  be 
genuinely  new  results.  They  are  useful  for  modeling  and  characterizing  acoustic  signals 
propagating  in  a  shallow  water  environment  characterized  by  high  attenuation  and 
transverse  isotropy.  This  is  an  environment  where  a  perturbation  treatment  of  the  bottom 
properties  applied  to  mode  summation  signal  synthesis  will  lead  to  erroneous  results.  The 
solution  of  equations  (35)  and  (36)  are  used  to  represent  the  anelastic  modes  in  terms  of 
the  elastic  modes,  permitting  a  detailed  analysis  of  the  physics  of  strongly  bottom 
interacting  acoustic  propagation.  The  effects  of  transverse  isotropy  and  attenuation, 
including  attenuation  induced  dispersion,  are  properly  accounted  for.  Stable  well- 
behaved  numerical  algorithms  exist  for  solving  the  complex  generalized  eigenvalue 
problem,  even  in  cases  where  the  the  matrices  involved  are  near  singular.  Our  next  step  is 
the  numerical  implementation  of  Eqs.  (35)  and  (39). 
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